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This report is submitted in accordance with the reporting provisions of 
the subject contract and describes the research effort and its accomplish- 
ments for the contract period. This study was performed by personnel of the 
Computational Mechanics Section of the Loclcheed-Huntsville Engineering Center 
(Lockheed-Huntsville) , Huntsville, Alabama, and Computational Mechanics Com- 
pany (COMCO), Austin, Texas. The Technical Representative of the NASA- 
Langley Contracting Officer was Dr. G.C. Olsen, Aero thermal Loads Branch, 

Mail Stop 395. 

1 . 1 BACKGROUND 

Finite element numerical methods are currently being used in computer 
codes for solving practical fluid flow problems. The advent of the super- 
computer is one of the primary reasons for the success thus far. However, 
future problems are destined to be more complex and will no doubt tax even 
the fastest machines. In conjunction with the release of the next generation 
of supercomputers (Cyber 2xx/GF10) , more powerful numerical algorithms will 
also be needed. Current methods utilize a grid of points to discretize the 
continuum which are fixed a priori and not changed during the computation. 

In addition, the order of the method, direction of differencing, and damping 
models, are all chosen by the code user. 

The success of finite element and finite difference codes often depends 
on the user's ability to discretize the domain and/or selectively increase 
the order of the finite element shape function to capture strong gradients 
within the domain. Currently, this requires an a priori knowledge of the 
location and strength of sharp gradients that occur in the flow field. Even 
then, obtaining optimal discretization and interpolation is a lengthy and 
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costly iterative procedure. Strong flowfield gradients occur in shock waves, 
expansion regions, and viscous Layers. The accurate determination of these 
regions is vitally important in determining the aerothermal loads on aero- 
space vehicles in supersonic flight. Body heating rates are particularly 
sensitive to the resolution of thermal gradients at the vehicle surface. 

The next generation of finite element methods to impact the computa- 
tional mechanics community will be the "self-adaptive" kind. In these ad- 
vanced methods, logic is built into the code to choose the grid of points, 
move them around, choose the degree of approximation, and generally adapt 
itself to the physics of the flow. Not only does this provide more reliable 
and accurate results, but it frees the (non-expert) user from making these 
decisions before running the code. 

1.2 OBJECTIVE 


The objective of this contract is to develop new computation methods 
for aerothermal heating analysis which utilize adaptive strategies. The new 
methods will be tested initially in trial codes and then implemented in 
Lockheed's GIM/PAGE code. Finally, a test problem will be run and compared 
with experimental data for code verification. 
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2 . WORK ACCOMPLISHED 


Adaptive procedures may be placed in one of three basic categories: 

1. Moving Meshes . The number of grid points is fixed, and the mesh is 
distorted so as to improve the quality of local approximations of 
the flow field and its gradient. 

2. Mesh Refinement (h-method) . The mesh is refined (i.e., the number 
of elements is increased, their dimension decreased) so that local 
accuracy is improved. 

3. Subspace Enrichment (p-method) . The local order of the approxima- 
tion is increased to provide a more accurate solution. In finite 
element methods, the mesh is fixed while the local degree p of the 
polynomial shape functions is increased. 

Regardless of which category an adaptive procedure falls into, it generally 
follows the steps shown in Fig. 2-1. The term structure refers to the basic 
mesh topology, the number of nodes and cells, the local order of the approxi- 
mation, the numerical scheme, etc. It is the framework within which the 
solution is obtained. Using an initial structure, a solution is computed. 

The "goodness" of this solution is then determined. A measure of "goodness" 
can be obtained by computing a posteriori error estimates. The measure of 
solution "goodness" can also include such things as the cost of the solution 
in dollars and the manhours required to obtain solution. If the "goodness" 
criteria is met then a solution of a specified "goodness" has been obtained. 
If the "goodness" criteria are not met, then the structure of the mathe- 
matical approximation is changed in some rational manner. This may involve 
moving nodes, adding more nodes and cells, and increasing the local order of 
the approximation. A "better" solution is now computed. This process is 
repeated until the "goodness" criteria are met. 
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This report consists of a collection of papers which document the 
various research efforts undertaken during this multiyear contract by 
Lockheed-Huntsville and COMCO personnel. 

Appendix A documents recent advances in error estimation and adaptive 
methods for finite element calculations. 

Appendix B documents the adaptive mesh strategy which is employed in 
several test codes as well as the GIM/PAGE code. 

Appendix C documents the implementation of a class of adaptive pro- 
cedures for time-dependent Euler equations in two dimensions. 

Appendix D documents implementation of an adaptive procedure which uses 
triangular elements and a FEM-FCT numerical scheme. 

The implementation of the GIM/PAGE code with adaptivity is documented 
in Appendix E. 

Development of a three-dimensional adaptive procedure is covered in 
Appendix F. 
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3. REMARKS AND CONCLUSIONS 


Adaptive finite element methods will have a significant impact on comp- 
utational fluid dynamics in the future. This report shows that adaptivity 
can be coupled to several numerical algorithms. Existing flow solvers can be 
enhanced with adaptivity. There is much work which needs to be done in the 
general area of overall adaptive strategy optimization. This involves the 
integration of both software and computer to realize an efficient analysis 
tool. Within a software/ computer structure, a particular adaptive strategy | 
may produce the least computationally expensive answers. This same adaptive 
strategy within another software/computer structure may perform very poorly. 
The type of data management technique may effect how an adaptive strategy 
performs within a particular software/ computer structure. In summary, 
additional work should be done to determine the effect of the software/ 
computer structure on an adaptive strategy. 
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RECENT ADVANCES IN ERROR ESTIMATION 
AND ADAPTIVE IMPROVEMENT OF 
FINITE ELEMENT CALCULATIONS 


J.T. Oden, P. Devloo,'and M.'Howe 


Texas Institute for Computational Mechanics, 
Department of Aerospace Engineering 
and Engineering Mechanics 
The University of Texas at Austin 


Abstract . We collect in this article a synopsis of methods and results’ on adaptive finite element 
methods. We outline methods for constructing a-posteriori error estimates for linear and nonlinear 
problems in mechanics. Adaptive methods are described and a variety of numerical results are 
given on applications to problems in fluid mechanics. 


1. INTRODUCTION 

How good are the answers? What can be done to improve them? These questions arise with 
increasing frequency among users of modern computational mechanics codes. They are 
fundamental, in that such questions relate to the basic goals of computational mechanics: the use of 
computational methods and devices to simulate mechanical phenomena. Yet much of contemporary 
research in computational mechanics is concerned with a myriad of other issues which, important 
as they may be, do not consciously and directly focus on those primitive and fundamental 
questions. When one does focus on those queries, a sequence of natural constraints are met that 
have a profound effect on the way one approaches the development of modem codes, numerical 
schemes, algorithms, and data management techniques for computational mechanics applications: 

Modulo natural deficiencies in the ability of the mathematical model itself to capture real 
physical behavior, we translate the first question into one that can be managed in mathematical 
terms: how accurate are the numerical solutions? The only plausible and general approach toward 
answering this question is to construct a-posteriori error estimates; i.e. to use the results of an 
initial calculation to estimate the local error in a finite element / finite difference approximation. 

Having obtained an indication of "how good the answers are," one can proceed to the second 
question: what can be done to improve them? The answer is clearly to use adaptivity of the 
approximation in some way: to change the structure of the approximation to improve accuracy, 
where by "structure" we mean the basic mesh topology, the number and location of nodes and 
cells, the local order of the approximations, etc. 

As is well known, there has emerged in the literature several methods for effectively :iltering 
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this structure: h-methods, in which the mesh is automatically refined to reduce error, p-methods, in 
which the local polynomial degree is increased; r-methods, in which a fixed number of nodal 
points are redistributed to reduce error; and combined methods, in which h - p, r - h, r - p,- 
combinations are employed. A survey of the recent literature on such adaptive methods has been 
compiled by Oden and Demkowicz [20]. 

What is especially significant about these answers to the basic questions , is that they have a 
great impact on the design of computational mechanics codes. To implement a rational adaptive 
scheme one must obey the following criteria in designing a programming strategy: 

1. Mesh Independence. Since the mesh itself may well be changing as the solution 
evolves, it is necessary to have Schemes which can be implemented on arbitr^ unstructured or 
quasi- structured meshes. This first criterion makes obsolete most existine body-fitted 
coordinate schemes common in finite difference literature. 

2. Robustness. Since the structure of the approximation is continually changing in an 
adaptive scheme, adaptive methods must be very stable under changes in mesh size, under mesh 
distortions, etc. 

3. Mathematical Basis. Since a -posteriori error estimates are necessary for an effective 
adaptive scheme, it is necessary that a solid mathematical basis exist for the adaptive methods. 

4. Geometry Independence. Modem computational methods, adaptive or not, must be 
able to cope with solution domains of arbitrary, complex geometry. The "real world" problems 
encounter^ in applications seldom have simple geometries for which many classical methods work 
well. 

5. Supercomputing. The significant data management problems inherent in adaptive 
strategies must lend themselves to supercomputing strategies-vectorization, parallelism, etc. 

6. Efficiency. Hopefully, when all features of an adaptive strategy are optimized in a 
program/computer structure, an efficient analysis tool will emerge. It is not necessary that the final 
product be capable of analyzing a given discretization as "fast" as possible; rather, the objective is 
overall optimization: to produce the best possible answers (in some sense) for a fixed level of 
computational effort. 

In our opinion, it is very clear that only finite element methodologies can fulfill all of these 
criteria. 

In this paper, we shall outline several recent advances in developments of the basis 
components of adaptive methods. We do not attempt to provide a thorough review of the literature, 
as this has already been the subject of a recent paper [20]. Rather, we provide summary comments 
in a few areas that we think stand out as important advances in the field. Naturally we are most 
familiar with our own efforts in this field, so we comment more fully on some of our own results. 

Following this Introduction, we give a brief summary of a few recent advances in adaptive 
finite elements. This is followed by several sections on general ideas behind a-posteriori error 
estimation, h-method data management, algorithms for fluid-mechanics applications, and some 
new results on numerical experiments with our adaptive codes. Finally, we comment on future 
directions of research in this field. 


2. RECENT .\DV.ANCES 

The state-of-the-an in adaptive finite element methods is adequately summarized in the volume 
of collected '-vorks and presentations m,ade at fhe Lisbon conference of 1985 These have recently 
appeared under the editorship of Babuska, Zienkiewicz, Gago. and Oliveira [ij. Here one v,ill fine 
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information of the basic ideas of h, p, r methods together with numerous applications to problems 
in solid mechanics and fluid mechanics. 

More recently a number of signigicant advances have been made which should be brought to 
view in the area of elliptic problems, we mention the important theoretical work by Guo and 
Babuska [14] on h-p methods. It is known that one can generally achieve a faster increase in local 
accuracy using p methods than h methods. By this it is meant that greater accuracy can be achieved 
with fewer degrees of freedom by increasing the local order p of the polynomial than by refining 
the mesh. This does not necessarily mean that the p methods offer a superior approach to solving 
elliptic problems, for one must add to this equation the significant factor of a data management 
scheme, which is often the life and death of an adaptive method. Babuska and his co-workers have 
shown, however, that the best possible approach to the accuracy problem, one leading to 
exponential convergence, is to simultaneously refine both h and p. The h-p methods have shown, 
in certain example problems, to produce exceptionally accurate results. At this writing, most of 
these results have been confined to one-dimensional problems and to linear elliptic problems in 
two-dimensions. There would appear to be some computational difficulties in extending these 
methods to time dependent problems, since there one must cope with the difficult issue of 
consistent mass matrices, stability and space-time approximation. However, it is possible that these 
difficulties may also be overcome with additional research. 

A production finite element code based on p methods is now being promoted and sold. This 
is the PROBE code, and its successful implementation of the p method has already an impact on the 
design of linearly elastic structures, see [25]. The simple r methods produced by Diaz and 
Kikuchi, and Taylor [12] have been used effectively in classes of problems in which one wants to 
keep the number of degrees of freedom more or less constant. In particular, in problems such as 
metal forming simulations, where one must solve a large number of nonlinear partial differential 
equations, it is natural to try to achieve the best possible accuracy for a fixed number of nodal 
points. Some simple moving mesh algorithms have been proposed which are easy to implement 
and which apparendy work well in two and three-dimensional problems. These have proved to be 
very effective for nonlineaf problems in plasticity in nonlinear solid mechanics. 

In general, moving mesh methods suffer from one defect: for a fixed number of nodes and 
fixed degree polynomial within each element, there is an inherent threshold of error which cannot 
be eliminated. Thus, with the exception of the work of Miller on moving finite element methods 
and the work mentioned above by Diaz and Kikuchi on r methods, most of the recent work on 
adaptive methods has focused on h-methods and p-methods. 

Perhaps the most significant recent advances in adaptive finite element methods have come in 
the area of time dependent problems. We mention in this regard the important work of Flaherty and 
his co-workers (see, for example, [8]) who have developed effective numerical methods for certain 
classes of parabolic problems. Additional references on this subject can be found in these papers. 
We also mention the construction of adaptive characteristic Petrov-Galerkin methods by 
Demkowicz and Oden [9,10] which involves not only the construction of the local a-posteriori 
estimates but also the construction of near optimal schemes for nonlinear convection diffusion 
problems with small diffusion coefficients. These results have recently been extended to solve 
Euler equations in two-dimensions. [27] . 

One area in which adaptive methods appear to be making some in-roads is in supersonic gas 
dynamics and general fluid mechanics. Several effective numerical schemes have been proposed 
by Lohner, Morgan and Zienkiewicz [16,17,18], and the authors [19, 21]. These schemes have 
been used effectively to solve two-dimensional steady state and transient problems in compressible 
fluid mechanics. 

More recently', Oden, Strouboulis and Devloo [19,23,24] have extended these methods to 
fluid mechanics in which moving domains are encountered. In particular, adaptive schemes have 
been developed for classes of problems in which flow interaction occurs due to the motion of one 
body or another through a flow field. Initial results on the application of adaptive methods to 
supersonic rotor-stator problems have produced some impressive results, some of which are 
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outlined later in this paper. These include examples in which adaptivity has resulted in a mesh with 
nearly 70 percent fewer elements than the uniform fine mesh required to produce equivalent 
accuracy. 

We comment on some of the components of an effective adaptive scheme below. 


3. A-POSTERIORI ERROR R.STIMATION 

The great majority of results on a-posteriori error estimation that have appeared in the literature 
in recent years is restricted to linear elliptic problems; however, a great deal of precision and depth 
of results has been possible for problems in this class. In [21.24,9], we have described a general 
method for a-posteriori error estimation that is applicable to broad classes of linear and nonlinetir 
problems, including parabolic and hyperbolic problems. Successful use of our method in 
determining error estimates for the finite element approximation of the Navier-Stokes equations has 
also been made [24]. 

■\n outline of the general method is provided by the abstract linear problem: 

Find u e V such that a (u,v) = f(v) V v s V (3.1; 

where a ( ■ , • ) is a bilinear form on V x V, V being a Banach space, and f is a linear functional on 
V. The Galerkin approximation of (3.1) in a finite-dimensional subspace of V is characterized 
by the problem. 


Findu^^ e such that a(u^, v^) = f(v*^'i V v*’ e (3.2'i 

We suppose that V — > H = H* V*, the inclusions being dense and continuous, for a 
Hilben (pivot) space H, V* being the dual of V, etc. If < . , . > denotes duality pairing on V* x V, 
then we generally have 

a(u, v) = <Au, v> 

If V = [v e V ! A v e H], then we also have <Au,v> = (A u,v)p| ; u,v e V where C , •) 
is the inner product on H and .\ e £ { W ^ , H.) 

In general, for finite element approximations, the form a ( • . • ) can be expressed as the sum 
of contributions from an assembly of E subdomains: 

u.v € V: a(u,v) = Z ((Au,v) -t- T (u.vi] 

e 

where ( • . ‘ ) denotes the H- inner product defined on restrictions of u and v to subdomain 
(element) e and T^(u.v) is the bilinear concomitant associated with boundarv’ terms on the boundar\' 
of subdomain e. • 

Let e^ = u - u^ denote the error and suppose that 


II v il-^ a(v.v) 


A-4 



sup a(w,v) 


llwll^ = 


Then 


!le^ 


V e V II V 11^ 

sup ale^'.v) 

V e V II V 11^ 


= II V 11^^ { S (r*’, v)^ + Fg ( , v) } 

V e V e 


(3.3) 


where r*' = Au - Au** = f - Au*’ is the local residual. To eliminate r^, we construct a local auxiliars- 
problem, for a function 0 defined by 


a( e «. V ) = R ( V ) ; e=l,2,...,E 


(3.4) 


where (v) = (r**, v)^ + F^Ce*’*, v) e*** being some appropriate approximation of e*' on the 
boundary. Setting Al^ the restriction of A over Qg and 

110*^ ll“A.c = 3e (0®. 0*") = < A'e 

introduce (3.4) into (3.3) to arrive at the a-posteriori estimate, 


2 1/2 

lle"llA SlZlie'll^ J 


(3.5) 


The functions 0^ are local error indicators. Of course we do not wish to solve the E equations 

to obtain the 0^. We are, thus, content to construct an approximate solution to (3.4) over some 

enriched subclass of functions so as to produce approximations 0^*^ of 0^. Several different 
methods of a-posteriori error estimation may result from different schemes for approximating (3.4). 
Alternatively, if one can derive local a-priori bounds such as 110® 11^^ < C II R® ll*,then (3.5) can be 
rewritten in terms of the residual functional R®. 

In many nonlinear problems, a step such as (3.3) may not hold, and instead, we bound the 
residual. For example 

,, u„ sup <Au-Au^v> sup fV/ h ^ T' / h M 

llr"IU= llvll {Z(r", v) + F (e", v)} 

vgVIIvII veV 


< { ZII0I|2 }l/2 


A,e 


We conclude this section with several remarks. 


L These examples provide global a-posteriori error (or residual) bounds in terms of local 
;rror indicators. By a special construction of test functions, truly local error estimates can be 
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obtained. For example, Demkowicz and Oden [9] studied a special Petrov-Galerkin method for the 
problem - e u" + u = f, and showed that the local error must satisfy the sharp a-posteriori estimate 




h2 


h- + en- 


|l r*> II 


l2(Q^) 


where is, again, the element residual. 

2. For a time -dependent problem, such as 
Jfj ( Ou/^t + A(u)) V dx dy = Jq f V dx dy 

tor arbitnu 7 test functions and v, and linear A, tlie fact that the error must be the function 
= u • leads, by direct substitution, to the evolution equation, 

Jq ( 3e*’/9t + A(e^)) v dx dy = - Jq r*’ v dx dy 

Thus, using a higher order approximation E** of e^ than that used in approximating u^, we arrive 
naturally at a system of equations for the evolution of error, 

ME-i-KE = R (3.6j 

Various dynamic error estimators can be constructed depending on how one constructs the 
approximation of e*' . In (3.6), M is the usual mass matrix associated with the approximation 

E^ = Ej Ej (t) Yj (x), E is the vector of nodal errors E-i, K is the stiffness matrix, and R the residual 
vector. 

3. For certain classes of problems, it is possible (or, at least, it may be assumed to be 
p<)ssiblej to obtain an estimate. 


Vv'-eV^I' 

•vhere t^e special class of local test fuentions used in approximating the local auxilliary 

problems (?.4). Then II 9^ - II Ac ^ estimated using standard results from finite 

element interpolation theory (see Oden and Carey [22]). In particular, if 9^ is the interpolant of 9^ 
over obtained using polynomials of degree < k, for an n - dimensional problem with 
quasi-uniform mesh refinements. 


0 .. - 0^ I 


m. q. 


< C h 


n/q - n/p -r k -f- I - m 


9.1 

^ k +■ l.p. n. 


(3.7) 


w ith • n, j Q the ^ ( O ) - seminorm 0 < p < ’o. and q - p / (p - 1 For the case m = 0, k = 
1 , p = q = 2, we obtain 


C',n 


< C h„- i 0.. ; 


L- (il^) 


2 ^ 


♦ .V. .V.i .11 — — V, U — i. 


^ 1 ^ _ 


'A C fliiVc* 
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iee-ehidfl-h^^i0e-eh 


AVG 


ch^^ie^ 


l,«o,e 


These estimates can judge the quality of the approximations of the local indicators, provided a 
means for computing estimates of the seminorms I * k+1 p c developed. 


4. FEATURES OF AN ALGORITHM APPROPRIATE FOR ADAPTIVE FEM 

Earlier in this paper we lifted criteria for the development of adaptive finite element codes for 
complex problems in solid and fluid mechanics. In this section, we summarize features of an 
adaptive code we have developed for two dimensional problems in compressible gas dynamics in 
which we have attempted to meet most of these criteria. 

4,1 Preliminaries 

We consider the motion of a perfect gas flowing through a two-dimensional domain Q cr ir 2_ 
If U = U (x,t) is the vector of conservation variables with p the mass density, m the linear 
momentum and e the total energy, it satisfies the following weak initial -boundary value problem: 

Find U e V such that 


+ Q (U) : V <|) ) dQdt -I- /iIq ()>(•, 0)dQ = J J F'^pdsdt (4.1) 

0 Q 3t n 0 dci 

V(j)€ W 

Here Q (U) is the Euler flux tensor. 



mj 1 

m 2 


p'* m^^ p(U) 1 

p‘l mjm2 

Q(U) = 

p-» mim2 1 

p-lmj2 +p(U) 


_p'l mj (e p (U) 1 

p*l m 2 ( e + p (U)) 


p(U) = (Y- l)(e-p-l m-m/2) 

where p is the thermodynamic pressure and y is the ratio of specific heats. 

Moreover, 

V = { V = { Vj, V 2 , V 3 , V 4 )'^ I Vi = Vi (x. 0 € L“ (0, T; L* (Q)); i = 1,2,3, 4) (4.3) 

W = { w = { Wj, W 2 , W 3 , W 4 }'!' I wj 6 C* [Q, T]), w • (x, T) = 0 ; i = 1,2,3,4} (4.4) 


F is the actual prescribed flux through the boundary 3Q and the following notation is used 
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T i 

U -^ = Z U- — 

0=1 3t 

Q:V4>= Z I Q^. _ 


i=l 0=1 


9X; 


Let us now consider an arbitnu7 time interval [x j, x 2I C ( 0, T] and modify the space of test 
funcrions to include functions which do not vanish at the final time, namely; 


W 




= { w = { wj, W2, W3, W4)’^ I Wj e C^Q X [ X j, X 2 1 ) ; i = 1, 2,3,4) 


Then we can state the weak-statement of the conservation laws over the space time subdomain Q \ 
1 X |, X 2 1 as follows: 

Find U e 2 such that 

In (U'r(-.t2)()> (.,X2) dQ= Jq (U'T(-,Xi)(1) (-,Xi) dQ. 

+ Q : V 4) ) dQ dt J F'^'bdydt V0e WM''^2 (4.5) 

Tj £2 9t 

Here V i’ 2 is appropriately defined as the solution space over the strip Q x [ x j, x 2 ]• 


I 


4.2 Solution Algorithm 

We obtain a finite element approximation of (4.1) by partitioning the space-time domain Q x 
fO.T] into subdomain Qx[t„, tj,^^] (with 0 = t q < tj < ... < tj, < t < ... t;,^- = T) by 
discretizing each subdomain and by employing (4..^ using the discrete spaces of test and tnal 
functions defined by the discretization. Moreover, by approximating the space-time integrals using 
numerical integration we get the following scheme [19): 

1. First Step: 

For each element Q^., compute such that. 


^OgdQ = 1 q U"^dQ - At/2 If^div Q (U") dQ 


(4.6) 


II. Second Step: 

Calculate such that. 



dQ = + At Q ( IT-’ " 


) : T Ph dQ 
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(4.7) 


-AtJ (t>7(Q(Uh'^l)-Q(U")) 

an “ 


an 


Here we assumed d^Jdi = 0 (i.e. the spatial grid remains fixed), we let F = Q n and denote 

Equations (4.6), (4.7) define a two-step TG/FELW (Taylor-Galerkin/Finite-Element-Lax 
WendrofO methc^ which has been introducted by Donea [13], studied by Baker et al. [21, refined 
by Lohner et al. [17] and others ( [3], [19]). The second step of the scheme, as given in (4.7) 
involves a global calculation of the form: 

M ( U = ( R ) (4.9) 

Here M denotes the consistent mass matrix, (R) the load vector whose definition can be easily 
deduced from (4.7) and (U) = [ Uj, U2, ..., U^]^ is the global vector of nodal 

unknowns. The inversion of the mass matrix can be performed by a Jacobi iteration [17] or a 
preconditioned Jacobi Conjugate Gradient [19]. 

The TG/FE-LW method provides us with a fast, multi-dimensional time stepping algorithm 
with a high resolution (high order of accuracy) in smooth regions of flow and which applies to 
unstructured adaptive grids. It is well known [17] that the algorithm suffers from a phenomenon 
of non-linear instability. To overcome this deficiency, artificial diffusion is added to stabilize the 
scheme in the presence of discontinuities ( [18], [ 19]). 


4,3 Flux-Corrected Transport 

The theory of Flux Corrected Transport has been developed by Boris, Book and others ( [4], 
[5], [6]) and it involves an attempt to systematically correct finite - difference transport schemes in 
order to avoid non-physical oscillations in the solution. Fully multi-dimensional FCT schemes have 
been constructed by ^esak [26]. Recently Lohner et al. [18] presented a flux-correction procedure 
of the TG/FELW scheme for systems of conservation laws. In this section we give a short 
exposition of the FCT - TG/FE-LW algorithm which we employed in some of our adaptive 
calculations. 

The FCT procedure consists of solving equation (4.9) by using a diffusion and an 
antidiffusion step. In the diffusion step a "strong" diffusion term is added to obtain a "transported 
and diffused" solution which is free of non-physical oscillations. In the antidiffusion step part, a 
"limited" amount of diffusion is subtracted from the right hand side (4.9) in order to steepen the 
solution at discontinuities and increase the accuracy in "smooth" regions of flow. 

In particular, we have: 

Step I: "Diffusion" Step 

Compute { from 

=(R) +{ V) (4.10) 

Here [V} denotes the vector of added diffusion with nodal contributions of the form; 
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(4.11) 


a<h,iT3U" 34),iT3U" 

''i = Io®x— — + Dy 

3x 3x ’^3y 3y 

For a mesh of quadrilaterals we let 

Dx =Dy =cAg 

where c is a constant and denotes the area of element Qg. 

Step II : "Antidiffusion " Step. 


Compute {U"'*’' ) as the limit of the sequence of iterates { > = 1,2,3,... defined by: 


M, = I ( F) 

^ [i+1] td [1] 

F = ( Mr - M) - V 

^ [i] 


(4.12) 


Here Ml denotes the lumped mass matrix and I denotes the flux limiting function which may be 
defined appropriately in order to prevent oscillations in the solution. In our applications we used 
the strategy of Zalesak [26] and Lohneret al. [18] to compute 1 (F^jj). 


4,4 An h Refinement / Unrefinement Strategy for Steady- State Solutions of 
High-Speed Compressible Flow 

An adaptive procedure for steady-state solutions of equations of compressible gas dynamics 
involves the following steps: 

For a given domain a coarse finite element mesh is defined which contains only a number of 
elements sufficient to model the basic geometric features of the flow domain (see Figure la). Each 
element in the initial mesh is assigned a "level" equal to zero. Then a finer mesh is generated by a 
bisection process, indicated in Figure lb, in order to obtain an initial grid with the "group" 
structure. Note that when an element is refined a group of 4 elements is defined and each of the 4 
new elements has a level one unit higher than the "parent" element. 

1 . For a given finite element grid determine the steady-state solution. 

6 over all M elements in the grid. Let 
e ° 


3. We scan groups of 4 elements and compute 

^GROUP- 


2. Compute error indicators 
®MAX “ 

l<e<M 
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where n\ is the k-th element in group m. 

4. Error tolerances are given by two real numbers, 0 < a, p < 1. 


If 0,> p 0 

C 


MAX 


we refine element Qg by bisecting it into four new elements. 


®GROUP ^ “ ®MAX 

we unrefine group m by replacing this group with a single new element with the nodes 
coincident with the comer nodes of the group. 

5. Go to step 1. 


4.5 Numerical Examples 

In this section we present examples of adaptive calculations of steady-state solutions of 
problems in high speed compressible flow. The error indicator employed in the numerical 
examples is given by the normdized gradient of the density; 


0 . 



^Ph 

max 

i=l,2 9xj 


Ph 


where denotes an average value of the density of element Qg . 


4.5. 1 Supersonic Flow Over a 20° Ramp 

We consider the problem of a Mach 3 flow (with y = 1.40) over a 20° ramp. The gas enters 
with uniform flow conditions through the left boundary of the domain and develops an oblique 
shock at the root of the ramp. 

A coarse initial mesh with the computed pressure contours are illustrated in Fig. 2. Adaptive 
mesh results are shown in Figures 3 and 4 with one and two levels of refinement respectively. The 

constants for the adaptive scheme were chosen a = 0.05, P = 0.15. The FCT version of the 
time-stepping algorithm was employed with c = 0.125. The results compare well with the exact 
solution except for some small disturbances downstream which are due to the artificial stagnation 
point at the tip of the comer. A three-dimensional view of the pressure is shown in Figure 5. 


4.5.2 Supersonic Flow in Expansion Corner. 

In this example, the steady supersonic flow through a 10° expansion is studied. The inflow 
Mach number was selected = 6 with y = 1.38. Figures 6 through 8 show the meshes 
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Figure 2. Supersonic flow over a 20® ramp. 

Initial mesh and pressure contours. 
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Figure 3. Supersonic flow over a 20° ramp. 

Mesh and pressure contours obtained with one level of rcfinemenL 
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Figure 4. Supersonic flow over a 20° ramp. 

Mesh and pressure contoun obtained with two levels of refinemenL 
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Figure 5. Suoersonic flow over a 20° ramp. . . • u 

Three-dimensional 'dew of the converged pressure funcnon obtained with two 

levels of refinement 
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F:gure 6. Supersonic expansion around a 10° comer. 
Initial mesh and density contours. 
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employed in the calculation with the corresponding density contours. The results were obtained 

with the FCT scheme with a = 0.05, P = 0.15 and c = 0.125. Striking improvement in the 
solution is seen to result from the refinement procedure. 

-i.5.3 Asynchronous Time-Stepping Procedures 

In the algorithms described in the previous paragraph the global timestep At is determined as 
the minimum allowable time step in thfe grid, namely: 

cTac 

At = min At^; Atg= (4.14) 

e=l,...,M lul + c 

Here C denotes the C.F.L. number, c is the local speed of sound in the element and luf = u"i + u^2- 

From the definition (4. 14) we see that since At ~ C ~ hg the timestep may be governed 

by the smallest element in the mesh. This choice of At guarantees stability and time-accuracy of 
the scheme. For steady-state calcumations however time-accuracy is not important and it may be 
more economical to employ asynchronous time-stepping by prescribing local time-steps. 

Let us denote by At j"°4e ^he nodal timestep of node j which is computed by the minimum of 

the time-steps of the elements which are connected to node j. Then, an Asychronous TG/FE-LW 
scheme may be employed as follows: 

I. First Step: 

For each element Qg, compute Ug"+^ such that: 

Ug"+’/2 = Uh"da.Ate Jo div Q(Uh")dD (4.15) 

c c 'y c 


IT. Second Step: 

Calculate such that, 

h j=l - ‘ 



+ At",o4e Q ( : V (()j dQ ) 


- At J (t)y(Q(Uu"+l/2).Q(u"))ndY 

+ At"o^®J Q ( Uh") n dy j = l,2, .... N (4.16) 

J 9QJ 
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We now demonstrate some of the features of the asynchronous time-stepping scheme using 
two numerical examples: 


4.5.3. 1 ■ The Reflecting Shock Problem; 

The statement of the problem is given in [19]. Figures 9 and 10 show the steady-state density 
contours obtained with the time-accurate and asynchronous algorithms, respectively. We note that 
the steady state was obtained after 130 time steps with the time-accurate scheme and after only 100 
time-steps with the asynchronous scheme, which represents 30% of savings in computational 
effort. 

4.5.3. 2. NACA 0012 Airfoil in Supersonic Wind Tunnel: 

We also considered the problem of a NACA 0012 airfoil in a supersonic wind tunnel with 
inflow Mach number = 3, y= 1-4 [19]. Figure 11 presents a comparison between the 
steady-state density contours obtained with the two schemes. The time-accurate scheme requires 
585 time-steps to converge while the asynchronous scheme converged after 496 time-steps. 


4,5.4. Elevon Cove Problem 

The Elevon Cove problem has to do with supersonic flow past a complex swan-like geometry 
of a portion of the space shuttle elevon. The problem is described in [3]. Figures 12 and 13 show a 
preliminary calculation of the problem with our adaptive Euler code. The mesh shown does not 
correspond dto a later unrefined mesh. This mesh is not yet optimal, since the program was still 
attempting to compute a new mesh at the time calculations were stoppeid. 


5. Features for an Adaptive Finite Element Algorithm for Transient Calculations 

We now present an example of an h-refinement / unrefinement strategy for transient 
calculations. The basic steps of the algorithm arc: 

a) Advance the solution N time steps. 

b) Do the following until no more elements can be refined: 

( 1 ) Compute the element error indicators 6g. 

(2) Refine all elements with 0g > Qjviax 

(3) Integrate the last N time steps with the updated (refined) mesh 

(4) Go to (1). 

c) Compute the element error indicators 0^ and unrefine all groups with 

®‘"gROUP - “ ®MAX 

d) Go to a). 

We note that the "do loop" in step b) converges when no more elements can be refined (the 
maximum level of refinement is fixed). Although the iteration in step b) guarantees a "fully 
updated" mesh it may lead to an expensive scheme if more than a few passes are required for 
convergence of the "do loop". A cheaper alternative is presented by the following "two-pass ' 
scheme: 

a) .Advance the solution N time steps. 
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b) Compute the clement error indicators 0^. 

c) Refine all elements with 0^ ^ P 

d) Integrate the last N time steps with me refined mesh obtained in c) 

e) Compute the element error indicators 0g and 

1 ) Unrefine all groups with 

^'"grouP - “ ®MAX 

2) Refine all elements with 

n Go to step a). 

In the following, we present two examples of adaptive refinement for transient problems. 


5 ! Rotating Cone Problem fill 

We consider the following advection problem; 


^ + div (aU) = 0 


j 0, r > 150 

U(x,y,0) = l ^ 

I 250[1 +COS — ],r< 150 

Here, r^ = + (y - 250)^ is given by the vector a (R, 0) = (R cos 0, - R sin 0) where R, 0 

are the polar coordinates indicated in Figure 14. 

This problem has been solved by many authors and it is considered as a benchmaric problem 
for algorithms for advection problems ([15], [11]). Here we show some results obtained with an 
adaptive SUPG algorithm [11]. Figure 15 shows some "fully updated" meshes which are obtained 
with the scheme outlined in the beginning of this section. For more details the reader should refer 
to [11]. 


5.2 A Problem of Supersonic Rotor-Stator Interaction 

We applied the "two-pass" adaptive algorithm to a problem of supersonic rotor-stator 
interaction. We consider now two rows of doubly-parabolic airfoils with thickness to lenth ratio 
equal to 0.08. Figure 16 shows some of these airfoils and the initial finite element discretization of 
the domain. We assume that the stator and the rotor have the same number of airfoils and we 
perform the computation on domains corresponding in one rotor and one stator airfoil while the 
presence of the remaining airfoils is simulated by periodic boundary conditions. In the figures the 
domain of the rotor airfoil is drawn twice. 

In the Figs. 17 through 25 we give the results of a supersonic calculation obtained with a 
dynamically adapted grid. The distance between consecutive airfoils of the rotor (and stator) is 
assumed equal to twice the airfoil length while the distance between the tail of the stator and the 
front tip of the rotor airfoil is taken equal to 0.2 of the airfoil length. We impose boundary 
conditions of supersonic inflow on the left boundary of the stator with the dependent variables 
equai to 


p = 1.4, pu = 4.2. pv = 0. pe = 8.8. 
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The inflow boundary conditions correspond to a free stream Mach number equal to three. 
Boundary conditions of supersonic outflow were assumed on the right boundary of the domain of 
the rotor. The steady-state solution which is obtained by keeping the airfoild fixed was used as 
initial condition. 

We have chosen P = 0.19, a = 0.06 and we defined the group error indicator to be equal with 
the maximum element error indicator of the elements in the group. We did not specify N but instead 
we revised the mesh every time the fifth nodes [23,24] of the rotor mesh coincided with comer 
nodes of the stator mesh (this resulted in mesh revisions every 10-12 time steps). We also note that 
all interface elements have been refined beforehand with the maximum level of refinement to 
facilitate the application of the sliding interface algorithm. 

In order to capture shocks of variable strength we used the "normalized" error indicators given 
in (4.13). It becomes clear from the numerical results that variable shocks are captured well and the 
mesh evolves dynamically to adapt to the solution of the rotor-stator problem. 

Results are shown in Figures 17 - 25. The initial mesh is that shown in Fig. 16. The first 
adaptive calculation for a steady-state initial condition is shown in Fig. 17a. The corresponding 
computed pressure profiles are shown in Fig. 17b. Note the symmetry of the shock lines, the 
continuous pressure fields across the mesh interface, and the fact that both unrefinement and 
refinement of the mesh were required to achieve the accuracy limits specified. The rotor blades are 
then allowed to move with unit speed and the mesh is dynamically refined. Plots are shown of 
calculated adaptive meshes and pressure profiles for 1/8, 2/8, 3/8. 4/8, 5/8, 6/8, 7/8 and 1 cycle 
('period) of the motion, during which a rotor blade makes a complete revolution from its initial 
position in Fig. 17 back to the same position. 

Several features of the computed meshes and solutions are noteworthy. In the initial 
steady-state case, only 16 unrefined elements appear. The size of these large elements, indicating 
small local error, is limited in the present caluclations by the distance from the tips of the rotor and 
stator blades: two elements in the present case since there must exist a sliding interface between 
them. A minor program modification could allow much larger elements in regions of small error. 
For the transient case, the number of larger elements (indicating substantial unrefinement) 
increases, and these regions of low error migrate over the mesh as solution evolves in time. 
Conversely, substantial refinement of the mesh is indicated at the interface and along shock lines. 
The method successfully captures shock interactions and the increasing density of pressure profiles 
downstream from the moving blades. The ratio of the number of elements in the adaptive mesh to 
that in the uniform fine mesh varies in time, but is typically 4,(XX) / 12,5(X), a reduction of 68 
percent! The initial coarse mesh of Fig. 16 contains around 4000 cells and is incapable of delivering 
the required accuracy, a fact not easily realized without an expensive computation. 


6. Future Directions 

We believe adaptive finite element methods will have a significant impact on computational 
fluid dynamics and computational structural mechanics in the future. These techniques, together 
with the modem parallel and array processors, will make obsolete many of the more popular 
methods in numerical analysis in use today. In particular, use of the body-fitted coordinate 
techniques, splitting methods such as ADI, etc. will probably lose some of their popularity since 
they arc not well suited for problems with unstructured meshes. 

It is likely that large gains are to be made in three-dimensional problems. Here more than 
anywhere else, one needs to do computations on a near optimal mesh where only a minimum 
number of degrees-of-freedom is required to produce a given level of accuracy. 

It is likely that new advances in parallel and array processing will bring the p-methods and h-p 
methods to the forefront, since, at least from a theoretical point of view, array processors may have 
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j Figure 10. Reflecting shock problem. 

j Adaptive grid and density contours obtained with the asynchronous time-stepping 

scheme. 

J 
J 


A-25 





Figure 1 1. NACA 0012 airfoil in supersonic wind tunnel. 

(a) Steady-state density contours obtained with the time-accurate scheme. 

(b) Steady-state density contours obtained with the asynchronous time-stepping 
scheme. 
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Figure 12. Eleven Cove Problem. 

Adaptive Finite Element Grid. 
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Figure 17. 



Supersonic flow interaction between rotor and stator airfoils: 

(a) Initial adaptively refined mesh for steady-flow through rotor-stator configuration. 

(b) Pressure contours for steady-flow through rotor-stator configuration. 
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Figure 19. Supenonic flow interaction between rotor and stator airfoUs. 

(a) Adaptively refined mesh at 2/8 of the rotor cycle. 

(b) Pressure contours at 2/8 of the rotor cycle. 
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Figure 20. Supersonic flow interaction between rotor and stator airfoils. 
( 3 ) Adaptively refined mesh at 3/8 of the rotor cycle. 

(b) Pressure contoun at 3/8 of the rotor cycle. 
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Figure 21. Supersonic flow interaction between rotor and stator airfoils. 

(a) Adaptively refined mesh at 4/8 of the rotor cycle. 

(b) Pressure contours at 4/8 of the rotor cycle. 
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Figure 22. Supersonic flow interaction between rotor and stator airfoils. 

(a) Adaptively refuted mesh at 5/8 of the rotor cycle. 

(b) Pressure contours at 5/8 of the rotor cycle. 
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Figure 23. Supersonic flow interaction between rotor and stator airfoils. 

(a) Adaptively refined mesh at 6/8 of the rotor cycle. 

(b) Pressure contours at 6/8 of the rotor cycle. 


A-37 



Figure 24. Supersonic flow interaction between rotor and stator airfoils. 

(a) Adaptively refined mesh at 7/8 of the rotor cycle. 

(b) Pressure contours at 7/8 of the rotor cycle. 
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Figure 25. 



Supersonic flow interaction between rotor and stator airfoils. 

(a) Adaptively refined mesh after one complete rotor cycle. 

(b) Pressure contours after one complete rotor cycle. 
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the capability of handling signiflcant local data management problems that might be associated with 
an implementation of p-methods. It is conceivable that a clever use of p-method philosophy in this 
computing environment may prove to be very effective for a wide class of problems. 

There is another point of view that is emerging from the adaptive literature: that is that two 
general types of a-posteriori estimation can be used in effective adaptive procedures. In one case, 
only a rather crude error estimator may be satisfactory to establish trends irt mesh adaptation that 
will lead to improved accuracies. The effectivity indices for such methods may not be close to 
unity, so that the actual error predicted may be quite far from the true error that exists in the 
approximate solution. Nevertheless, a scheme may result which is truly adaptive, in the sense that 
the actual local error is systematically reduced below some threshold. Parallel to these methods are 
methods in which a great deal of sophistication is used in an a-posteriori error estimation. Here, 
with additional expense, quite accurate estimates of local errors can be obtained. This leads one to 
speculate that there may emerge in the future adaptive codes with two or three levels of 
sophistication: one in which an adaptive scheme is used to produce a near optimal solution for a 
fixed level of effort; secondly, a post-processing operation in which very precise estimates of the 
local error are produced and presented, perhaps in terms of error contours, for residual evaluation 
toward obtaining a final evaluation of the quality of the solution. Again, if this quality is not 
acceptable to the analyst, he may choose to re-run the problem through additional adaptive cycles to 
produce furhter improvement in local solution quality. 
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Appendix B 


The adaptive mesh strategy to be described is an h-method applied to 
hexahedral elements for three-dimensional and quadrilaterals for two- 
dimensionals wherein the mesh is refined or unrefined (coarsened) when a 
solution quality test function falls outside preassigned upper and lower 
bounds. For clarity the two-dimensional strategy is described. It differs 
from the three-dimensional strategy only in the number of elements which 
comprise a group (4 vs 8) and the number of new sub-elements created during 
a refinement (4 vs 8). A set of "adaptation" rules are listed which are 
used to implement this strategy. 


1.1 General Description . The adaptive mesh strategy involves the following 
steps. 

1. For a given domain L, such as that shown in Fig. B-la, a coarse 
finite element mesh and an initial solution are available. 


2. As the adaptive process will be designed to handle groups of four 
elements at a time, a finer starting grid is generated by a 
bisection process, indicated in Fig. B-lb, to obtain an initial set 
of element groups. 

3. The adaptive procedure is initiated by computing solution quality 
indicators r^ over all M elements in the grid. Let 


’^MAX' 


,max„ 
1 e H 


r 

e 


4. Next, scan groups of a fixed number P of elements and compute 


k _ P 
GROUP" k=l 
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Fig. B-1 (a) A Coarse Initial Mesh Consisting of Four-Element 

Groups and (b) The Refinement and Unrefinement of a 
Group of Elements 


B-2 


where e is the element number for group k, P = 4 for two- 

IV 

dimensional grids and P = 8 for three-dimensional grids. 

5^ The solution quality bounds are defined by two real numbers, 0 a, B 
1. If 


r 

e 


< B r. 


MAX 


element r is refined. This is done by bisectinR r into four 
e e 

new sub-elements. If 

k 

r ar 

GROUP = MAX 


group k is unrefined by replacing this group with a single new element 
with nodes coincident with the comer nodes of the group. This is 
always possible because each group is itself the result of an initial 
bisectioning. 

This general process can be followed for any choice of a solution quality 
indicator. 

2.2 Data Structures . An important consideration in all adaptive schemes is 
the data structure and associated algorithms needed to handle the changing 
number of elements, their node locations and numbers, and the element labels. 

As noted in the preceding paragraphs, the algorithm is designed to pro- 
cess (refine or unrefine) in groups of four elements at each local refinement/ 
unrefinement step. Consider, for example, the case of an initial mesh of 20 
square elements shown in Fig. B-2. Assign to each element in this mesh an 
element number, NEL = 1,2, . . , ,NELEM and to each global node a label NODE. The 
array, NODES(J,NEL) relates the local node number J(J = 1,2, 3, 4) of element 

NEL to the global node number NODES. In addition, the coordinates X ,Y 

j «J 


B-3 










of each node are also provided relative to a fixed global coordinate system. 
File these numbers in two arrays, 

NODES(J,NEL) = the array of global node number 

assigned to node J of element NEL 

XDO(JCO,NODE) = the array of JCO — coordinates of 
global node NODE (JCO = 1 or 2). 

If a solution quality indicator signals that an element should be refined, 
say element 11 in the example, some system for assigning appropriate labels to 
the new elements and nodes must be devised. Toward this end, a convention can 
be established that defines the connectivity of the specified element with its 
neighbors in the mesh. This information is provided by a third connectivity 
array. 


NELCON(NC,NEL) = the connection of element NEL, 

where NC = 1 , 2 , . . . , 8 

As seen in Fig. B-2, each side of an element may be connected to two other 
elements so that NELCON is dimensioned thusly; 

NELC0N(8,MAXEL) 

with MAXEL an appropriately large number. 

The entire refinement process (or its inverse — the unrefinement process) 
just described is accomplished by specifying a series of element levels. For 
example, the initial coarse mesh could be assigned level 0. When an element 
is refined, its sub-elements belong to a higher level, level 1, and when these 
sub-elements are refined, elements of level 2 result, and so on. In this way, 
if the maximum level any element in the mesh can achieve is limited, then the 
maximum number of elements the mesh can contain is also limited. In general, 
no such limit need be set. 
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Thus, the bookkeeping of element and node numbers evolved in a refinement 
process is monitored by the arrays NODES(.,.), XC0(.,.), NBLCON(.,.), and an 
array LEVEL(NEL) which assigns a level number to element NEL. Initially, the 
same level can be assigned to all elements, and this level is an arbitrary 
parameter prescribed in advance by the user. Thus, provisions are now in hand 
for an arbitrary, dynamic renumbering of elements and nodes. 


2.3 Adaptation Rules . Several rules must be established to successfully 
implement the refinement or coarsening of a mesh. The following "element" 
rules are employed: 


1. An element may be refined only if its neighbors are at the same 
refinement level or higher. 

2. If a "neighbor" element of an element to be refined is at a lower 
level of refinement, it must be refined first. 

3. Refinement of an element results in creation of eight sub-elements 
for three-dimensional and four sub-elements for two-dimensional 
meshes . 

4. To be eligible for coarsening a group of elements must not contain 
another group of elements and each element of the group to be 
coarsened must not be connected to a "neighbor" element of a higher 
level. 

For example, if element 11 if Fig. B-2 is to be refined, we proceed through 
the following steps: 

1. An intermediated node is common to two members of a group only. 

2. An intermediate node that is created along a domain boundary cannot 
be constrained. 

3. If an element and its neighbor both of which are at the same level 
are connected to a third element at a lower level, then the 
intermediate node which exists along the edge common with the third 
element is constrained. 

4. If a group of elements is eligible for coarsening, then the 
intermediate constrained node along the edge common to an element 
which is not a member of the group will be eliminated. 
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5. If a group of elements is eligible for coarsening, then the node 
along the edge common to this group and its neighbor group will 
become constrained. 

6. If a group of elements is eligible for coarsening, then the 
intermediate node along a domain boundary edge is eliminated. 

Use of the above rules can be illustrated by considering the uniform grid 
of four elements shown in Fig. B-3a. Suppose element A is marked for refine- 
ment. By applying element rules 1 and 3, element A is divided into sub- 
elements, I, II, III, IV as shown. Application of node rules 1 and 2 dictates 
that the nodes marked by circles be constrained. Nodes marked X are un- 
constrained. 

Next, let element III be chosen for further refinement. Element III 
cannot be refined since one of its neighbors, B is at a lower level. Refine- 
ment of element III before element B would violate element rule 1. Therefore, 
element B is refined as shown in Fig. B-3b. Note that node B is no longer 
constrained, since node rule 2 no longer is satisfied. Node Cl remains 
constrained. 

Now that element B has been divided into elements V, VI, VII, VIII, 
element rule 1 can be applied. Figure B-3c illustrates this division. 

Suppose the group of elements V, VI, VII, VIII shown in Fig. B-3c is 
marked for coarsening. This group is not eligible for coarsening until the 
group of elements, a, B, Y, w has been coarsened. Element VII has neighbors B 
and w which are a higher level. This violates element rule 4. 

Now, let the group of elements, a, B, Y, w be marked for coarsening. 
Element rule 4 is satisfied and elements a, B, Y, w are replaced by element 
III. The intermediate constrained nodes associated with elements a, B, Y, w 
are eliminated through use of node rule 4. The intermediate node along the 
upper domain boundary is eliminated using node rule 6 . 
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I, Fast Reflnement/Unref inement and Moving 
Mesh Methods for Unstructured Meshes 


J.T. Oden, P. Devloo, and M. Howe 

Texas Institute for Computational Mechanics, 
Department of Aerospace Engineering 
and Engineering Mechanics 
The University of Texas at Austin 


ABSTRACT . New adaptive finite element methods are presented for the 
analysis of unsteady invlscld compressible flow in arbitrary two-dimen- 
sional domains. The procedures described herein are used in conjunction 
with a semi-explicit two-step algorithm for solving the time-dependent 
Euler equations in two space dimensions. Two schemes are presented for 
monitoring the evolution of error, and error estimates are used as a 
basis for a mesh refinement strategy. The capability of unrefinement 
(adaptively coarsening the mesh) is also included. The methods do not 
require a structured mesh and are*. applicable to quite general geometries. 
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1. INTRODUCTION 


Many would agree that the most fundamental and important questions 
facing users of modem computational methods for flow predictions are the 
following: 

I. How good are the answers? 

II. How can one obtain the best possible answers for a fixed 
computational effort (or a fixed computing budget, fixed 
manpower level, or a fixed and limited computing 
capability)? 

The first question is exceedingly difficult since it Includes both the 
Issue of the validity of the physical and mathematical model of the flow 
phenomena Itself as well as the Issue of the quality of the numerical 
approximation of the equations characterizing the model. To simplify 
matters for purposes of the present discussion, we shall dispense with the 
first Issue and take for granted that the classical Navler-Stokes or, in 
the present paper, the Euler equations are adequate models of nature for 
the applications In mind. Thus, the first question reduces to a word; 
accuracy — how accurate are the numerical solutions? 

The second question Is seldom asked, but it is Intrinsically con- 
nected to the first. It is common practice in applications of computa- 
tional fluid dyanmlcs to the complex flow domains, to generate extreme- 
ly fine finite difference meshes In hopes of capturing all important 
features of the flow, even though the location of these special points 
of interest changes In time. This leads some to employ fine meshes In 
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all positions of the flow domain where some Important aspect of the 
flow might possibly manifest itself. The quality of results is gen- 
erally judged by the invariance of solutions to further refinement if 
one can afford the cost of another calculation. The fact that coarse 
mesh solutions may be adequate In much of the domain at most Instants 
of time cannot be exploited In traditional fixed mesh schemes. 

After some thought about these Issues, rather broad answers to 
the fundamental questions present themselves: 

I. Accuracy . To determine the accuracy of a computed 
solution, one can attempt to develop reliable a-posterlorl estimates 
of local error. In other words, one might hope to be able to develop 
procedures which use the evolving computed solution to determine sharp 
estimates of local -errors In various norms over each mesh cell and at 
each time step. 

II. ’’Optimal*' Meshes . Use adaptive procedures to contin- 
ually change the structure of the mesh — the size of mesh cells, the 
location of grid points — so as to keep the local errors within a 
preassigned limit. 

Obviously, the second answer assumes that one has some means to 
measure the local quality of the numerical solution and, therefore, 
presumes the availability of some-- type of a-pqsterior*i error estimate. 

We describe. In this paper, algorithms and results developed in 
an attempt to more sharply resolve these answers, particularly that to 
question II, for a class of problem In compressible flow. More speci- 
fically, we describe here a class of very effective adaptive schemes 
for time-dependent Euler equations In two dimensions which employ both 
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mesh refinement (when the local error Is large) and mesh "unrefinement" 
(when the local error is small) and which generate the appropriate mesh 
changes as the solution evolves in time. This requires that we estimate 
the local approximation errors at each time step. However, only an Indl- 
catlon of the relative error between successive meshes Is essential In our 
methods; the issue of very sharp a-po8teviori estimates of local error (ou 
answer to question I) Is one of great concern to us and Is the subject of 
other papers [8,18,19]. 

In designing an adaptive scheme for Euler equations, we keep the 
following guidelines In mind: 

(1) Unstructured Grids . The method must be virtually grid 
Independent and global-coordinate free. While an initial 
mesh can be defined to model the basic geometry of the flow 
domain and the initial data, thereafter it must be possible 
to automatically add or eliminate cells and grid points as 
needed to monitor local accuracy levels. This requirement 
considerably lessens the attractiveness of body-fitted 
coordinates, many elliptic/algebraic mesh generators, and 
various factorization algorithms which exploit such regular 
mesh topologies. 

(2) General Geometries and Boundary Conditions . The method must 
be applicable to arbitrary flow domains with virtually 
arbitrary geometry, general in-flow and out-flow conditions, 
and general boundary conditions. 

(3) Solid Mathematical Basis . Since, by Its nature, any sound 
adaptive method must employ some type of local error 
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estimator. It Is important that the methods employed have a 
reasonably firm mathematical basis, e.g., that a-prior*i or 
a-posteviori error estimates exist and that the convergence 
characteristics of the method are acceptable. 

(4) High Accuracy . The method should be capable of delivering 
high-order accuracy. 

(5) Robustness . The method must be numerically stable and not 
sensitive to singularities, distortions in the mesh, or 
Irregularities in the data. 

(6) Super computing . The method should lend Itself to modem 
supercomputing methods for accelerating computational speed, 
such as easy vectorlzatlon or implementation on parallel 
processors, etc. 

(7) Computational Efficiency . The method and the supporting 
algorithms and data structures must be computationally 
efficient. 

We feel that these criteria can be best met by finite element methods. 
In the present work, we use as the basis of our adaptive schemes a semi- 
explicit method used by several other authors (e.g., [19,4,15,16,21]): a 

two-step Lax-Wendroff /Taylor-Galerkin scheme. It is far from optimal (and 
does not satisfy all of our critCTla) , but is perfectly adequate to use in 
conjunction with our adaptive scheme. Schemes which fulfill all of these 
criteria are under development and will be reported in subsequent papers. 

We remark that there is a growing literature on adaptive methods in 
computational fluid mechanics. Adaptive procedures for incompressible 
viscous flow problems were developed by the authors in a series of recent 
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papers (see, e.g., [8,18,19]). These methods employed a variety of differ- 
ent adaptive strategies, but did not come as close to satisfying the above 
criteria as the methods discussed in the present work. The general subject 
of adaptive finite element methods is dealt with in a forthcoming volume of 
collected papers edited by Babuska, Zlenkievlcz, Gago, and Oliveira [2]. 

For a survey of adaptive finite difference schemes, see the works of ' 
Anderson [1]. Also, Berger and Ollger [4] and Berger and Jameson [5] have 
recently developed adaptive finite difference methods for hyperbolic 
conservation laws. Still other types of adaptive methods for hyperbolic 
problems have been recently proposed by Demkowlcz and Oden [10,11]. 

Following this Introduction, we develop weak formulations of a class 
of problems in compressible gas dynamics. These space-time formulations 
are shown to be the. basis of a class of Lax-Wendroff /Taylor Galerkin 
schemes. Our derivation of this family of algorithms is nonstandard, in 
that we show that a two-step scheme follows easily from the use of a 
numerical quadrature scheme for evaluating appropriate flux integrals. In 
Section 3, finite element nodels of the space-time formulation are 
introduced, and in Section 4, we discuss the important issue of a 
posteriori error estimation. Section 5 of this paper is devoted to a 
detailed discussion of adaptive strategies. These include an h-method, 
wherein the mesh is refined or unrefined when local errors fall outside a 
preassigned upper and lower bound, and an r-method, in which the mesh is 
automatically distorted to equldlstrlbute error. In Section 6 of the 
paper, we present the results of several numerical experiments on 
two-dimensional problems. These results illustrate that the performances 
of the adaptive schemes are quite acceptable for a class of complex flow 
problems. 
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2. PRELIMINARIES 


We consider the motion of a perfect gas flowing through a domain n 

over a time interval [0,T1. We shall confine our attention to two- 

2 

dimensional cases* Q cz K ; we denote by D the space-time domain* 

D = n X (0,T) and by the boundary of fi . The motion of the gas Is 

governed by the global balance laws of physics and the second law of ther- 
modynamics. Thus, If U = U(x,t) , (x,t) 6 D , Is the 4-vector of conser- 
vatlon variables, U = {p , m , E} , with p the mass density, m the 
linear momentum, and E the total energy, and If d^ and dS denote 
Lebesque measures of area (volume) and length (area) of Q and 3Q 
respectively, then we demand that U satisfy the following system of 
conservation laws:-. 



( 2 . 1 ) 


Here, Q(U) is the flux and n Is the unit outward normal to 312 . If 
m^ , m^ denote Cartesian components of m , then 

T 

U = {p , iDj , n>2 , E} 


Q(U) = 


p~^ + p(u) 

-1 

p 

p"^mjE + p(U)) 


* p m.in^ 

I -12 

j p m~ + p(U) 

I p"^(e + p(U)) 


( 2 . 2 ) 


n = {nj^ , n2)^ ; P(U) = (y “ 0(E - p ^ n • m/2) 
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In (2.2), p Is the thermodynamic pressure and y is the ratio of specific 
heats, assumed here to be constant. In addition to (2.1), U must satisfy 
the entropy production Inequality 

~ I pn(U) dn + I n • (mn(U) + e“^ gjdS & 0 (2.3) 

■'fi ~ ■'an “ “ ~ 

with n(U) the entropy density of the gas, 0 the absolute temperature, 
and q the heat flux, as well as an Initial condition, 

U(x,0) = Uq(x) , X 6 n (2.4) 

where Is given. 

It Is of fundamental Importance to note the smoothness requirements on 
U In order that (2.1) make sense mathematically. Conservation laws (2.1) 
hold when the components of U are bounded measurable (with respect to 
Lebesque measure In x ) functions on D . Thus, we may seek solutions in 
the function space 

V = {V = {Vi, V^, Vy = V^(x,t) 

g l“(o,T ; L^(J2)] ; 1 = 1, 2, 3, 4} (2.5) 

In particular, (2.1) Is not equivalent to the classical Euler equations, 

Uj. + dlv Q(U) =0 (2.6) 

3Q ,/3x.) since solutions of (2.1) may 
al 1 

not possess derivatives across surfaces In D , However, the conservation 
laws and Initial conditions are fully equivalent to the following weaTc 
boundary-initial value problem: 

Find U G V such that 


(with = 3U/3t and dlv Q = J 
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(2.7) 


f + 9<y> • 

Jn *' 


uj $(‘,0)df2 = A J dS dt 

n ^0‘'3n~ 

for all i 6 W 

where F Is the actual prescribed flux through 30 and W Is a suitable 
space of test functions; e.g.> 

W » <t>2* *^’3» I '^1 “ <J»^(x,t), 

6 C^(D) , ^^(x,T) = 0 . 1 = 1, 2, 3. 4} (2.8) 


In (2.7), we use the notation 

d» - y u ; 

- a 3t ’ 

05*1 


2 4 34) 

9 : ’ i - I I Q„i ^ 

i=*l a==l i 


On the other hand, if U Is known to have Integrable derivatives in 
D everywhere except on a family of surfaces > then we may 

consider the problem 

Find U G V such that 

f (?t ^ dt - I y^(-,o) 6(*,o)do 


- I f {Sfc lull - h n!l}dS dt 

k=i Jr, 
k 

+ I uj ^ dO + lo L ^ ds. dt 


-H 

^0 ha 


F ^ dS dt for all J 6 W 


(2.10) 



Here, are the speeds of propagation of discontinuities across and 

V and W are appropriately redefined, e.g,, 

V = {V = (Vj, V^, V^, V^) I 6 L*(D) , 

Vi S - D - r^j) 

W = {$ I -l>i 6 C°(D) , (t^(x,T) = 0} 

where Fj^^^ are the surfaces on which suffers a jump. If F is not 

a prescribed flux but Is merely a notation for Qn , then these flux terms 

cancel and do not appear In the formulation. 

Consider an arbitrary time Interval [ 12 * 12 ] d [0,T] and Include 

in W functions c^(x,t 2 ) * 0 . Let to be a subset of n such that 

ti) o u F. = 0 > and let F = Qn . Then another weak statement of the system 
^ Vi - . 

conservation laws over to x [t 2 »T 21 iss 
Find U G such that 

(-U^ + (div Q)^ j)dn dt 

(U 

f if " $(-,T2))dn -0 (2.11) 

J 0 ) 

for all ^ G W“’"^ 

with and appropriate spaces of trial and test functions. 

Remark . It Is well known that .(2.7), (2.10) and (2.11) may all' 
possess non— physical solutions since none of these formulations Involve the 
entropy Inequality (2.3). Thus, In general, we seek solutions to (2.7), 
(2.10) or (2.11) in the subset K c: V ; K = {v G V , v satisfies (2.3) for 
appropriate 0 , q(V0)} . 
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3. FINITE ELEMENT APPROXIMATIONS 


Finite element approximations of the gas dynamics problem are obtained 
by a direct approximation of (2.10) or (2.11) on finite-dimensional spaces 

I 

approximating the spaces V and W . The spatial domain 0 is partitioned 
into a collection of finite elements 0^ over which the components of 

trial functions V are approximated by polynomials of degree k . In this 
way, we construct a family of finite dimensional spaces of the type 

= {y^ * {Vj, V^, V^, 6 V I 

6 (0^) , i = 1, 2, 3, 4} (3.1) 

where P. (^^ ) is the space of polynomials of degree k defined over 0 

iC 6 • ' G 

Alternatively, we can use 6 » where space of 

e 

tensor products of polynomials of degree k on (e.g., is 

spanned by bilinear functions, Q2^^e^ biquadratics, etc,) . In addi- 
tion, a family finite dimensional spaces of test functions is 

also constructed. We then consider Galerkin approximations of (2.7), 

(2.10) or (2.11) by seeking solutions to these equalities in , with V 
and W replaced by and , respectively. 

3.1 A Two-Step, Lax-Wendroff/Taylor-Galerkln Scheme . We next derive 
a special semi-discrete, weak formulation from (2.11) which provides the 
basis for the construction of a popular family of finite element schemes. 

We proceed with the following steps: 
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i) Partition the time interval [0,T] according to 0 ■ t^ < t^ 
< t^ < . . . < = T ; 

ii) Apply the weak balance law (2.11) to a typical time Interval 
>‘n’ ■ 'n 

ill) Set " 0 in (2.11) suggesting the ultimate use of a 
time- invariant grid (we relax this assumption later) ; 
iv) Replace the time integrations in (2.11) by the elementary 
midpoint quadrature rule 


f n+1 
■'t 


f(t)dt ~ At f 


n+Jj 


At » t - t , = f(t • + At/2) 

n+ In n 


Thus, with 0 ) ® n , we obtain the semidiscrete approximation 
I ih dn = j y” dn + At I dn 


- At I 

Jan 

for all 


(3.2) 


where u” = U, (x, t ) , etc., U, being the approximation of 
"'ll n "*n 

U , and 9^^** is the flux at the half step, 

3"*'' - 8(3r‘') (3-3) 


v) To obtain an approximation Uj^ ' , we use (2.11) again for time 
Interval [t^, » this time replacing the time integrals by 

a simple strip rule and Integrating by parts the divergence terms 
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‘ n 

e 

^ e 



-4^1 J(dlv Q")d!) 

e 



for all 


thus 

arrive at the algorithm, 


1) 

With (u” , g'' = known at 

compute using (3.4) 

the n th time step. 

2) 

Compute using (3.3) 


3) 

Compute using (3.4) 


4) 

Go to 1) 



This algorithm is the finite-element based two-step Lax-Wendroff /Taylor 
Galerkin scheme (see [20,7]). It is one of a family of methods advanced by 
Donea [12], studied by Baker and Kim [3], and successfully refined and used 

by LOhner et al. [15,16]) and Bey et al. [6] in finite-element applications 

! 

in fluid dynamics. This semi-explicit method is of second order in time 
and can experience spurious oscillations near shocks and other types of 
irregularities in the solution. These deficiencies must be reckoned with 
in implementing the method, 

3.2 Artificial Viscosity . As noted earlier, artificial viscosity 
terms are usually added to schemes such as that employed here so as to 
dampen out oscillations in the numerical solutions near shocks. The 
calculations described subsequently were performed adding a Lapidus- 
viscosity term, which, at time step t^ , is of the form 

- 7 • (c(u) . (3.5) 
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where 


3u 


w v» , 

c.(u) “ C I -7—^ I (no sura on 1 ) 
1 ** dX j 


tll 

C is the Lapldus constant, = m^/p is the i component of the flow 
velocity, and a is a parameter which determines if viscosity is to be 
Included implicitly (a “ 1) or explicitly (ct = 0) . The viscosity term, 
written out in component form, is 


- 1 


h- O 5 e = 1, 2, 3, A 




Setting o “ 1 , we obtain for step 2 of the procedure [instead of 

(3.2)), 

4 -2 


Us- 


dn + At 


I I c.M u"*; ♦'> , 
a-i 1=1 ‘ “-i 


dn 


- At f jJ[c(u") • 7 u")n dS 

Jsn 

- 9"*'= = ! Ih 

- At I n)dn 


for all admissible test functions 


ih 


3.3 Details of the Finite Element Algorithm . The details of the 
implementation of the algorithm described above are crucial to successful 
computations. In this work, we use meshes of four-node quadrilateral (Q^ 
elements over which the components =,1,2, 3, 4) of U are piecewise 

bilinear functions. Similar approximations and algorithms are used by Bey 
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et al. [6]. In addition, so-called group approximations of the flux 


0^^(a = 1,2,3,^; 1 = 1,2) are employed so that these com.ponents are also 
piecewise bilinear functions determined by their values at element nodes. 
In general, this finite element approximation will be of the form, 

N 

j = l 

(3.7) 

"al ' 

J=1 


where N denotes the total number of nodes in the discretization, and 
, 0^ . are values of at node j , and . are the global 

Ct Ot 1. »» J 

piecewise bilinear basis functions. 

As noted earlier, we advance the solution in time in two steps. It is 
Important to note that the first step is essentially local, computed over 
each element, while the second is global and contains the artificial 
viscosity terms: 

First Step : For each element , calculate a constant element vector 

^ f rom 
a,e 



(3.8) 


Second Step : For each node j , calculate by solving the 

following system of equations 
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^6^ 



uJ.n+1 

a 


N 

I ( 

j = l J 


dn) 


j.n 



- At 


I,„ "8 ‘’"s *1 


ds 


C3.9) 


xn 


Here, Q denotes the elementwise averaged value of the flux. The 


coefficients t. are defined to be constant over each element, 

p 


T« “ cA 
B e 


3u 


3x, 


where c Is a global constant (c = 1 In the examples), A denotes the 

e 

area of , u^ denote the components of the fluid velocity. 

To speed up the calculation, we precalculate and store the following 
element Integrals before the time stepping Is started: 



dn , 


^i<tj dn , 



3x, 


dn 


3(J^ 3(t^ 

3^ 3Xg 


dn 


l,j = 1,2, 3, 4; 6 = 1,2, 3, 4 


An element-by-element Jacobi Conjugate Gradient method Is used to ob- 
tain the solution of the matrix problem in the second step. Due to the 
structure of the mass matrix, the iterative solver requires only a few 
iterations to converge fully. 
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3.4 Boundary Conditions . In the finite element schemes developed 
here, we implement the following three types of boundary conditions: 

(a) Supersonic Inflow . On the part of the boundary with 
supersonic Inflow, the values of all the conservation variables are 
imposed; 

(b) Supersonic Outflow . On the outflow part of the boundary, 
the values of the conservation variables and the normal flux are 
unknown. Boundary conditions of supersonic outflow are implemented by 
adding the contribution of the boundary integral of the normal flux to 
the right-hand side of the equations of the second step; and 

(c) Solid Boundaries . On a solid boundary, the normal component 

of the velocity u = u. n- is zero. We note that, in general, the 
n p p 

nodal directions are not uniquely defined. In such calculations, we 
compute the normal directions at the nodes which satisfy global mass 
conservation at the steady state, namely. 



1 r 

'I 

n n1 

B-1 



3.5 Hourglass Instabilities . We now show that the Taylor-Galerkin 
scheme presented above can propagate undetected spurious solutions. To 
demonstrate this, let us consider;. the scheme applied on the 2-D Burger's 
equations on a mesh of rectangular elements. Burger's equations may be 
obtained from the above formulation by redefining the flux as follows: 
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U..0, 



2 3 


3 

0 


0 


Consider a rectangular element with the following nodal solution at time 

t : 
n 

{u^}” = [0, 1. 1, 0]*= 

= to, 1. -1. O]*" 

{u^}" = (0, 1. 1, of 

“ to, 1, -1, 0]*^ 

Then the scheme gives, 

= to, 1, 0, 0]*" 

and, by letting c = 0 (no artificial viscosity), we get; 

= {u^l" 

This means that the scheme propagates "hourglass" solutions undetected. 
This fact explains why in the numerical examples the method produced oscil- 
lations near the outflow boundaries. This hourglasslng phenomenon can be 
eliminated by considering each quadrilateral as a patch of two triangular 
elements joined along one diagonal of a quadrilateral. 
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4. ERRORS 



The adaptive finite element methods described here Involve two basic 
components: 

1) Error estimation — the determination of a-posteviovi estimates 
of the evolution of error In the numerical solution; and 

2) Adaptation — the automatic restructuring of the approximation so 
as to reduce local element errors and the computational effort. 

In this paper, adaptive procedures are based on estimates of error In 
a single principal dependent variable, such as the density, pressure or the 
entropy. We shall choose the density p as the driving factor In adaptiv- 
ity, although other choices could be used In the algorithms developed here. 
Two basic procedures are used to estimate local element errors. 

4.1 Evolution Equation for Error . Consider the continuity equation 
for the evolution of mass density through a domain fi with known flow 
velocity u . A weak form of the continuity equation is 

p dfi = - V*(u p)(J) dn 

Q h 

for all 6 W (4.1) 

I 

A semi-discrete Galerkln approximation of (4,1) consists of seeking an 
approximate density such that, over some suitable finite-dimensional 

space of test functions W, , 

I \ pt " I 

for all (j>j^ 6 (4.2) 
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If W. C W , we may choose ♦ ■ In (4.1), subtract (4.2) and obtain the 
n n , 

h li 

following evolution equation for the error e (x,t) =» p(x,t) - p (x,t): 

I ell dO » I - 7‘(u dfi 

Jjl h t Jjj ~ h 

for all 6 (^.3) 

The exact and approximate solutions are related according to 

p = (4.4) 

where e^ Is the approximation error. Thus, the error satisfies the 
evolution equation, 

I ('^®g eN)dJl = 

for all ♦ 6 W (4.5) 

where <r,_,<t> is the residual functional, 
n 

<r, ,(^> = - I + V»u p^ij)) dfl (4,6) 

h Jjj ^ t ~ ~ 

If we replace by ® (4.3) and the evolution equation 

reduces to merely the orthogonality condition (4.3), which Is automatically 
satisfied by error. 

We obtain an approximate evolution equation for the error as follows: 
let denote a fine-grid approximation of e^; i.e., 

e^(x,t) ~ E^(x,t) = ^ E*^(t) 'l»«(x) (4.7) 

N 

where 'l'jj(j) denotes a polynomial basis function defined on a subgrid of 
finer mesh size than that used to calculate p^ . Then, Introduction of 
(4.7) into (4.5), and replacing (ji by gives 
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I ("nm e” * k(u)^ - 


rjj(t) - 0 


N = 1, 2, 


(4.8) 


where 


”nm = Jo '^N 


V’^M 


(4.9) 




Many possible ways for implementing (4.9) present themselves. These 
equations, for example, need not be global in the sense that an element-by- 
-element or patch of elements in a fine mesh obtained through a mesh re- 
finement may produce sufficient accuracy to allow for an adequate indication 
of the evolution of error. The local velocities u and residual r, can 

" u 

be Interpolated using Qj-approxlmations on a fine mesh level. Several of 
these alternatives are under study and are to be the subject of a forthcom- 
ing report. 


4.2 Interpolation Errors . Let u be a smooth function defined over a 
regular domain . The W^’^(il) -semi-norm of u is defined by 


w’^’P(n) i l+j=r '3xJ-3x^ 


i.j^O 


(4.10) 


where I ^ p £ " and r is a non-negative integer. 
The Sobolev norm of u is 
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lul 


W^»P(fl)u 


^ '"'"ip 

tk=o »r’P(n) 


1/ 


(A. 11) 


Let G be an arbitrary convex subdomain (a finite element) of G 

Mr 

over which u Is Interpolated by a function u^ which contains complete 
piecewise polynomials of degree k . Then, It can be shown (see Oden and 
Carey [17]) that the local Interpolation error In the Vp’^(G)-seml-norm Is 


n, 

” T.l“» 


w“’P(G) 


S C 


k+1 




n 

" w 




where 


(A. 12) 


h » the diameter of the domain G 

p s the diameter of the largest sphere that can be Inscribed Inside G 
n » the dimension of the domain G 
p' =» p/(p - 1) 

C 3 a constant Independent of h , P , and u . 

If p Is proportional to h and If It remains proportional In refinements 
of G defined by parametrically reducing h , we have 


n 

• 1 

- - + k+1 - 

m 


Im p 0 ^ 

m , p , (v 

P 


(A. 13) 

with 1* L _ p =■ 1* 1 _ _ 

mtP.G w“*P(G) 

t etc. and 

E = u - Uj^ . 

• 


Such estimates can be used to devise crude adaptive schemes. Suppose 
that u on the right side of (A. 13) Is replaced by a finite element 
approximation Uj^ and that 



Indicates that the local error in the W°*^(G) seminonn is proportional to 

. . ^ ,.n/p' - n/p + k+1 - tn • • 

the error indicator, h 

i) n=2, Tn*0, k=l, p = p' = 2 


ul. , . Some choices are: 

k+l ,p 


k*(G) ^ ^ 


In this case, one must approximate the ’*-semi-norm of u over G ; 
i.e,, the L^-norm of second partial derivatives of u . 


li) n = 2, p 


1, k = 0, m = 0 . 


® 'lI(G) = 


average ' 


“ Ch^ max|7*u(x) 
x6G 


Such estimates can give only rough indications of local errors in 
sufficiently fine meshes. However, they are usually easy to implement and 
our experience is that they can provide a very effective basis for mesh 
refinement strategies. 


5. ADAPTIVE MESH STRATEGIES 


Let us suppose that we can calculate an error indicator 0^ for each 
finite element in a given mesh at a time t . This indicator is, in 

general, a real number representing the local error in a suitable norm, and 
it is computed using one of the procedures described in the preceding sec- 
tion. The decision to adapt the numerical procedure (to refine the mesh or 
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to move nodal points) Is based on whether or not local error Indicators 
exceed preassigned tolerances. We shall describe two adaptive procedures 
In this section. 


5.1 An h-Reflnement /Unrefinement Method . Our h-procedure involves 
the following steps. 

1) For a given domain 0 , such as that shown in Figure la, a coarse 
finite element mesh Is constructed which contains only a number 
of elements sufficient to model basic geometrical features of the 
flow domain. 

2) As our adaptive process will be designed to handle groups of four 
elements at a time, we generate a finer starting grid by a 
bisection. process. Indicated In Figure lb, to obtain an Initial 
set of element groups. 

3) We Initiate the numerical solution procedures on this Initial 
coarse grid, and compute error indicators 0^ over all M 
elements in the grid. Let 


■^MAX 


max 0 
lSe<M * 


4) Next, we scan groups of a fixed number P of elements and 
compute 


0 . 


GROUP 


P 

I «. 

k»l k 


where ej^ Is the element number for group k . We take P = 4 
In our current codes. 
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5) Error tolerances are defined by two real numbers, 0 < a,B < 1 . 

This is done by bisecting 0^ into four 

''GROUP MAX 

we unreflne group k by replacing this group with a single new 
element with nodes coincident with the corner nodes of the group. 
This is always possible because each group is Itself the result 
of an initial bisectioning. 

This general process can be followed for any choice of an error 

i 

indicator. Moreover, it can also be Implemented at each time step in the 
numerical schemes discussed in Section 3. 

5.2 Data Structures . An Important consideration in all adaptive 
schemes is the data structure and associated algorithms needed to handle 
the changing number of elements, their node locations and numbers, and the 
element labels. 

As noted in the preceding paragraphs, the algorithm is designed to pro- 
cess (refine or unreflne) in groups of four elements at each local reflne- 
ment/unreflnement step. Consider, for example, the case of an initial mesh 
of 20 square elements shown in Figure 2. We assign to each element in this 
mesh an element number, NEL = 1 ,2, . . . ,NELEM and to each global node a label 
NODE. The array, NODES(J,NEL) relates the local node number J(J = 1,2, 
3,4) of element NEL to the global node number NODES . In addition, the 
coordinates X^,Y of each node are also provided relative to a fixed 
global coordinate system. We file these numbers in two arrays. 


If 


0 ^00 
e MAX 


we refine element 0 


new subelements. If 
^k 


O' 


£ Q0 
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NODES (J,NEL) 


the array of global node numbers 
assigned to node J of element NEL 
XC0(JC0,N0DE) = the array of JCO — coordlantes of 
global node NODE(JCO = 1 or 2) . 

Suppose that an error Indicator is computed that signals that an 
element should be refined, say element 11, In the example. We must have 
some system for assigning appropriate labels to the new elements and nodes. 
Toward this end, we can establish a convention that defines the connectivity 
of the specified element with Its neighbors In the mesh. This Information 
is provided by a third connectivity array, 

NELCON(NC,NEL) = the NC connection of element 

NEL NC = 1,2,..., 8 ; 

As seen in Figure 2, each side of an element may be connected to two other 
elements so that Dimension NELCON =* (8,MAXEL), (with MAXEL an appropriately 
large number) . 

The entire refinement (or its Inverse — the unrefinement process) 
just described Is accomplished by specifying a series of element levels. 

For example, the Initial coarse mesh could be assigned level 0 . When an 
element is refined, its subelements belong to a higher level, level 1, and 
when these sub-elements are refined, elements of level 2 result, and so on. 
In this way, if the maximum level any element in the mesh can achieve is 
limited, then the maximum number of elements the mesh can contain is also 
limited. In general, no such limit need be set. 

Thus, the bookkeeping of element and node numbers evolving in a re- 
finement process is monitored by the arrays NODES (•,•) , XCO(»,*) , 

NELCON(»,0 , and an array LEVEL(NEL) which assigns a level number to 
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element NEL . Initially, the same level can be assigned to all elements, 
and this level is arbitrary parameter prescribed in advance by the user. 
Thus, provisions are now in hand for an arbitrary, dynamic renumbering of 
elements and nodes. If, for example, for the mesh in Figure 2, if element 
11 is to be refined, we proceed through the following steps: 

(1) Loop over the neighbors of element 11 (which is 
made possible with the NELCON array) and check 
the level of the neighboring elements relative to 
the level of element 11; 

(2) If any neighboring element has a level lower than 
11, then the element cannot be refined at this 
stage; 

(3) If 11 can be refined (as is the case in Fig. 2), 
we generate new element numbers (thus changing 
NELEM and new node numbers for unconstrained 
nodes) ; 

(4) Compute the connectivity matrix NELCON for the new 
elements; 

(5) Adapt the connectivity matrices for the neighboring 
elements (since the refinement of 11 has now 
changed this connectivity) ; and 

(6) Interpolate the solution between the unconstrained 
nodes. 

It is clear that some strategy is needed to test if a designated element is 
appropriately connected for a refinement to take place. 
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Consider, for example, the uniform grid of four elements shown In 
Figure 3a and suppose that the error estimators dictate that element A Is 
to be refined. Thus, A Is divided Into four elements, I, II, III, IV, as 
shown, and the solution values at the junction nodes, shown circled In the 
figure, are constrained to coincide with the averaged values between those 
marked X . Note that the connectivities change In this process, e.g,, the 
connectivities 4 and 8 of element B are different. 

Next, assume that an additional refinement Is required, and that we 
must next refine element III. We Impose the restriction that each element 
side can have no more than two elements connected to It. Thus, before III 
can be refined, element B must first be refined, as indicated in Figure 
3b. The constrained node B1 in Figure 3a now becomes active, while node 
Cl remains a constrained node. With B bisected, we proceed to refine 
III into sub-elements a,B,Y»5» and new constrained nodes, again circled 
in Figure 3c, are produced. In this case, only element B had to be 
refined first in order to refine III, but, in general, the number of 
elements that must be refined in order to refine a particular element 
cannot be specified. The following subroutine determines the necessary 
refinements prerequisite to refining an element NELl: 

SUBROUTINE DIVIDE(NEL1 ,NEL2) 

NELl = the input element that needs to be refined 

NELl if NELl has 
been divided 

NEL2 = output element 

NELD = element that needs 
to be divided prior to NELl 
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Then, symbolically, we have the algorithm (for the example in Fig. 3), 


— Repeat 
NELl = III 

CALL DIVIDE (NELl, NEL2) 

— WHILE (NEL2.NE.NEL1) 
NELl = NEL2 

CALL DIVIDE (NEL1,NEL2) 

— END WHILE 

UNTIL (NEL2.EQ.III) 


5.3 Moving Mesh (Node Redistribution) Methods . Another family of 

adaptive schemes we have considered Is a node-redistribution scheme which 

progressively moves a fixed number of nodes as to reduce local error. One 

basis for such schemes is to equldlstrlbute error at each time step. 

For example, let 0 be an error Indicator for element fi In a mesh 

e e 

containing a fixed number M of elements in a two-dimensional mesh. Let 

h = h(Xj,X 2 ) be a mesh function such that 

h(x, ,x„) = h = dla(fi ) for (x, ,x_) G D 
1 2 e e i Z e 

and note that, approximately. 




(5,1) 


with = dXjdx^ (this being exact for domains which are unions of square 

elements). Let 0 = 0(Xj^,X2) be mesh function which gives the local error 

Indicator when evaluated at a point (0 = 0 for x G Q ), We wish to 

e - e 

minimize the total error indicator functional, 

Mr 

j( 0 ) = I 0 ^ dn (5.2) 

e=l h 
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subject to the constraint (5.1). Using Lagrange multipliers i this leads to 
the optimality condition, 


6(j + A(f h"^ dJ2 - M)) = 0 , 


or 


(s. H- • 0 

e 


or 


3 ^®e 


Suppose that meas(fJ ) = o h and that 0 is of the form 

e 0 e e 

0^ = h^f(u) . Then, Integrating this last result over a typical element 


gives 


'n 


oh^ 0. h®"^ f(u)dn = Xo h^ 

e e e o e 


Hence, the optimal mesh size distribution results when 


L 


0 dft =■ Xo /a = CONST, 
e o 


(5.3) 


In other words, 

.2 


Indicators 


0 


to obtain the optimal mesh, we must equidlstrlbute the 


J ® 

To use this result to redistribute nodes, we proceed as follows (cf. 
Diaz et al. [12]): -> 

1) Generate an initial (generally regular) mesh with a fixed number 
M of elements and compute a trial solution on this mesh at one time step; 

2) Compute the corresponding error indicators 0^ ; 
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Figure 4 , Calculation of area center-of-error 
N 

X to equidistribute element error 
indicators in a cluster of four elements. 
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3) For a group k of P elements (with P always 4 in this 

work), let A denote the area of element 1 In the group. The 
®i 

area-weighted Indicators for group k are the P-numbers, 



4) Let y denote a vector from the origin of a global coordinate 
®1 

system to the centroid of element e^ of group k . Then the center of 
error of group k is defined as the vector 



5) Relocate the node at the center of group k to lie at the vertex 

- k 
of X ; 

6) Continue this sequence of operations over each group h of four 
elements until the new location of each node does not change more than a 
preassigned tolerance. 

This process should approximately equldlstrlbute the element error 
indicators. 


6 . NUMERICAL EXAltPLES 

In this section, we present the results of several numerical 
experiments on representative test problems. Six examples are presented, 
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the first five involving steady-state solutions and the last a transient 
problem. In the steady-state examples, one of the following two strategies 
Is used: 

(A) A. 1. The numerical solution Is computed on a fixed mesh and is 

advanced in time until a steady state is reached. 

A. 2. After convergence to a steady state, error Indicators 0 

e 

are computed over each element. In the calculations dis- 
cussed below, we employ the interpolation estimates and use 

= A Jp^l^ ~ A I i9p^/3n|dS (6.1) 

e e ~ ® Jan 

e 

where A is the area of the element. 

A. 3. The mesh Is refined/unrefined using the criteria and 

algorithms discussed in the preceding section. 

(B) B.l. Same as step A. 1. 

B. 2. After convergence to a steady state, error Indicators 0^ 

are computed according to 


0 = A 

e e 


I. 


7p^ • 7p^ dfl 


( 6 . 2 ) 


B.3. In applying the node redistribution Cmovlng mesh) algorithm. 


a modified error indicator 0^ is employed which is designed 
to be always greater than unity even when ** ® par- 

ticular, we use 
a0 

®e “ ^ ^ Y0 


In our examples, a = 81 , B •» 1 , and y “ 8 . 
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B.4. Nodes are redistributed a total of K times using the 
procedure described in Section 5.3. In the examples, we 
take only two iterations (K » 2) . 

We proceed to the examples. 


6.1 Shock Reflection Problem . We begin with a problem for which an 
exact solution is known and which has been used as a benchmark problem by 
others. 

The problem involves the steady flow of a perfect gas in a rectangular 
duct in which density, velocity, and energy are prescribed in each of four 
triangular wedges in such a way that the appropriate jump conditions (the 
Rankine-Hugoniot conditions) are exactly satisfied. Thus, a problem of 
shock reflection for which an exact solution is known is obtained. Dimen- 
sions and data are given in Figure 5. In this and all the other problems, 
the solution is considered to have converged to steady state when the 
magnitude of the L^-norm of the density is reduced by three orders of 
magnitude. 

The time step is monitored by the formula 



2 yP I i2 2 *''2 

Here, C ^nd |u| ^ ^2 * Y ~ 1*^0 • The constants 

multiplying the artificial viscous terms were selected locally as: 


3u" 



3v" 

3x 

> 

' e 

T = A 

y e 

3y 


e 




Figure 5 . 


A shock reflection problem, 
variables are prescribed as 
outflow values are computed 
laws. 


Inflow values of the conservation 
indicated in regions I and II, and 
in III to satisfy the conservation 
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where the bar denotes average element values. A Lapldus constant of 1.0 
was used In all calculations. 

The results of a uniform coarse initial mesh approximation are shown 
in Figure 6. The computed density contours are also shown in this figure. 
Note that only a rough indication of the location of the shock is possible 
with this mesh. 

A much better resolution is given in Figure 7 where the adaptively 
refined mesh shown is computed with refinement parameters a = 0.10 , 

B = 0.50 (recall Section 5). Note that no "unrefinement” appears to have 
taken place with these parameter choices, but that the simple error 
estimation scheme is capable of detecting the general area of the shock 
line. The much improved density profiles are indicated in the figure. 

Still better results are obtained with the same a and 6 but with 
two levels of refinement, as Indicated in Figure 8. Note that in this case 
large elements appear in the mesh, indicating unrefinement as well as 
refinement, of the original mesh. The corresponding density surface is 
given in Figure 9 where quite sharp shock fronts are observed. Note some 
spurious oscillation is encountered near the outflow boundary, as should be 
expected from the deficiencies of the algorithm noted in Section 4. 

The same problem was also analyzed using the node redistribution 
algorithm discussed in Section 5.3 with 20 node redistribution iterations. 
Results are shown in Figure 10. There, the original coarse initial mesh of 
Figure 6 is progressively distorted to conform to the reflected shock 
locations. Corresponding density contours are also given in the figure. 
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l] 


Figure 7 . Reflecting shock problem. Mesh and density contours obtained 
with one Tcvel of refinement (a = 0.10, 3 = 0.50). 
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6.2 NACA 0012 Airfoil In Supersonic Wind Tunnel. In this example 


the supersonic flow through a narrow wind tunnel containing a NACA 0012 
airfoil is studied. The inflow Mach number was set at M,, =• 2 , with 
Y = 1.40 and symmetry is exploited to reduce the computational effort. 

The Initial coarse mesh and density computed contours are given In 
Figure 11, Note that the critical features of the solution — the 
reflected shock and contact discontinuity — are lost with this coarse 
mesh. A refined/unrefined mesh, obtained with parameters a « 0.10 , 

6 = 0.10 is shown in Figure 12 together with a greatly improved density 
approximation. In these and subsequent calculations, a CFL number of 0.5 
and a Lapldus constant of 1.0 were employed. Results of a node-redistri- 
bution scheme for the coarse mesh are shown in Figure 13. In these 
results, ten iterations of the node redistribution algorithm were used. 

6.3 Supersonic Flow in a Wind Tunnel with a Step . The steady-state 
solution of the problem of a wind tunnel with a s.tep introduced into the 
flow is next considered. The inflow Mach number was selected = 3.0 
and Y =* 1.40 • The initial coarse mesh is shown in Figure 14 with the 
corresponding density profiles, and results of the adaptive refinement/ 
unref inement scheme with ot = 0.15 and 6 = 0*20 are shown in Figure 15. 
The mesh refinement algorithm was also used, with the mesh and density 
profiles obtained after 10 iterations shown in Figure 16. We see that the 
adaptive scheme captures well the features of the flow including the con- 
tact discontinuity at the top near the point of reflection of the bow 
shock. However, some oscillations are present downstream, and they are 
believed to be due to the non-mcnctonlcity of the solution algorithm. The 
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results presented for the reflnement-unreflnement procedure have been con- 
strained by a maximum number of 2000 nodes or 2000 elements that can be 
allowed. In the refined mesh shown, this constraint has been achieved. 

6.4 Supersonic Flow Over a 20* Ramp . We next consider the steady 
supersonic flow through a conduit with a 20-degree ramp. The gas (with 
y - 1.4) enters as a uniform M = 3.0 flow through the left side of the 
ramp and a shock develops at the ramp root. A coarse initial mesh and the 
computed density contours are Illustrated in Figure 17. For this problem, 
a reasonably good indication of the orientation of the shock is obtained. 

Adaptive mesh results are shown in Figures 18 and 19 for choices of 
the parameters of a = 0.20 and 0 = 0.50- with one and two levels of 
refinement, respectively. Notice that spurious oscillations at the outflow 
boundary above the ramp root, due to the hourglass oscillations described 
in Section 3, cause unnecessary refinements in this region. Similarly, in 
regions between the shock and the ramps, some unnecessary refinement re- 
sults from oscillations in the numerical solution. Nevertheless, striking 
improvement in the quality of the solution is seen to result from the 
refinement procedure. 

In this particular problem, the node redistribution algorithm works 
remarkably well. A computed distorted coarse mesh, obtained after ten 
applications of the node redistribution algorithms, is shown in Figure 20 
with the resulting’ density contours. 

6.5 Blunt Leading Edge of 8* HTT Panel Holder in Hypersonic Flow . 

The problem of the blunt leading edge of the 8’ HTT panel holder in a 
supersonic flow field with freestream Mach number M^ = 6.57 , y = 1.38 
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Initial mesh and density contours. 


Figure 12 . NACA 0012 airfoil in supersonic wind tunnel 
Mesh and density contours obtained with one 
level of refinement (a = 0.10, B = 0.10) 
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Figure 13 . NACA 0012 airfoil in supersonic wind tunnel. 

' ^ Mesh and density contours obtained after 10 

’ j applications of the mesh redistribution algorithm. 

:j 
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j^igure 14. Supersonic flow in a wind tunnel with a 


step. Initial mesh and density contours. 




Figure 15 . Supersonic flow in a wind tunnel with step. 

Mesh and density contours obtained with one level 
of refinement (a = 0.15, 6 = 0.20). 
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applications of the mesh redistribution algorithm. 
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Figure 17 . Supersonic flow over a 20° ramp. 

Initial mesh and desnity contours. 
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Figure 19. 



Supersonic flow over a 20® ramp. Mesh and densitv 
obtained with two levels of refinement (a = 0.20, 


- contours 
g = 0.50). 
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and 0® angle of attack was solved to obtain the steady-state solution* 

This problem has also been studied by Bey et al. [6]. 

A coarse mesh solution is indicated in Figure 21 and an adaptively 
refined/unrefined mesh and solution, obtained for a = 0.05 and S » 0.15 , 
are shown in Figure 22. A distorted mesh and corresponding density map are 
Indicated In Figure 23. In this particular problem, neither the h-method 
nor the r-method gave particularly good results, as a poor approximation of 
the solution between the shock and blunt body results from spurious oscilla- 
tions in the basic time-marching algorithm. In the case of mesh adaptation 
using redistribution, the solution actually diverges after four passes 
through the adaptive scheme due to the badly graded (hourglassed) mesh 
produced from the oscillations of the adaptive scheme downstream of the 
shock. 


6*6 Transient Adaptive Solution for Supersonic Flow Over a 20^ Ramp . 
In all the examples presented above, a time-accurate time stepping scheme 
is used, but the adaptive scheme was not used stepwise for the transient 
solution since our primary interest was to increase accuracy in the steady- 
state solution* The adaptive method used to track transient fronts is 
described as follows: 

1) Choose a structured mesh with the finest mesh size to be allowed 
in the calculation to be the initial mesh. This is done to avoid 
large variations of the time-step during the time-stepping. 

2) Ever N time steps (N = 50 in the present problem) go through 
the ref inement-unref Inement process (only unrefinement after the 
first N time steps). 
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Figure 21 . 


Blunt leading edge in hypersonic flow field. 
Initial mesh and density contours. 
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Mesh and density contours obtained after, with 
one level of refinement ( = 0 . 05 , = 0 . 15 ). 
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Figure 23 , Blunt leading edge in 
hypersonic flow. Mesh 


and density contours 


obtained after 4 


applications of 




V/X 


the mesh 


redistribution 


algorithm. 
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The above adaptive strategy was employed to solve the problem of a 20* 
ramp which is suddenly introduced In a supersonic flow field with =» 3.0 , 
Y =1.40 . The solution was integrated to steady state, and it is demon- 
strated that the mesh adapts to the shock front as the shock front moves 
from its initial to its steady-state position. 

The initial coarse mesh is shown in Figure 24 and the evolution of a 
refined/unrefined mesh for various time intervals is illustrated in succes- 
sive figures. Figures 25-29. The refinement parameters used were o = 0.05 
and 8 = 0.25, and a total of 250 time steps were used to track the .solution 
from its initial to the final steady state. The final steady result is 
similar to that obtained earlier and shown in Figures 18 and 19. 



Figure 24 . Transient supersonic flow over a 20° ramp. Initial mesh. 
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Figure 26 . Transient supersonic flow over a 20“ ramp. Mesh and 
density contours after 100 time steps (a = 0.05, 

3 = 0.25). 



Figure 27. 



Transient supersonic flow over a 20° ramp. Mesh and density 

contours after 150 time steps (a = 0.05, 6 = 0.25). | 

L 
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Picture 28 . Transient supersonic flow over a 20° ramp. Mesh and density 
contours after 200 time steps (a = 0.05, 6 = 0.25). 


j 

J 
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Figure 29. Transient supersonic flow over a 20° ramp. Mesh and 
density contours after convergence to steady-state, 
(a = 0.05, 8 = 0.25). 
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Chapter 1 

This report concerns an adaptive finite element code, capable of solving 
transient and steady-state problems in compressible inviscid fluid flow. 
Unstructured triangular finite element meshes were used for the basis of the 
adaptivity. 

Some research has been carried out in the area of adaptivity (see [2],[4],[7] 
and [10]), although many different approaches have been pursued. For example, an 
approach involving quadrilateral elements was seen to provide accurate results in 
[4]. Devloo, [4], used one type of quadrilaterial element division, resulting in the 
generation of contrained nodes (i.e., the nodal solution is contrained in terms of 
other nearly nodal solutions.) This approach combined with a scheme (FEM-FCT, 
see Chapter 3), capable of capturing line discontinuities in the flow, was seen to 
improve the acciuracy of the results. Others, such as Bank [2], have opted for 
triangular elements and two different types of division. This choice of two divisions 
results in a mesh of unconstrained nodes, noted to provide accurate pictures of fluid- 
flow interaction. 

In all these references, adaptivity was seen to improve the solution and the 
speed with which the solution was computed. Considering the types of problems 
(i.e. compressible fluid flow), an adaptive triangular element approach with a 
suitable numerical scheme was hence adopted as the basis of this report. The 
adaptive part of the code forms the basis of the work carried out an an already coded 
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FCTG-FEM "solver" (that used by Strouboulis in [13]). Parts of the "solver" had to 
be changed to accommodate triangular elements. Full details of the adaptive part are 
presented in Chapter 4. 

Several classical fluid flow problems were analyzed. The initial conditions 
and solutions were obtained from references such as [11], [12], thereby enabling 
valid comparisons to be made. The problems all involved some form of line 
discontinuities such as shock/shock-interaction. The adaptivity was seen to help 
capture the true form of these discontinuities, providing similar solutions to those 
contained in [12]. The shock resolutions were accurate enough to justify the use of 
such a scheme. 

A brief summary of this report would thus be as follows: Chapters 2 and 3 
describe the FEM-FCT numerical scheme used; Chapter 4 describes the adaptive 
strategy; Chapter 5 includes numerical examples and; Chapter 6 gives suggestions 
and recommendations for further developments. 
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1 


CHAPTER2 

nNITE ELEMENT APPROXIMA’nONS 

Before describing Ae details of the Flux-Corrected Method (adapted to the 
Hnite Element Method - See Chapter 3), a brief background to the class of problems 
and method of solution will be given. 

All of the problems analysed involve the time-dependent flow of a 
compressible inviscid fluid in two-dimensions. If we denote the domain of flow by 
fi and the time interval of this flow by [0,T] then the true domain D is defined by 
D= ft X [0,T). A sketch of such a domain is given in Figure 2.1 including the 
boundary of the domain. This domain is then used to compute a solution vector U, 
consisting of the following four components : 


U = [ p , pu , pv , pe ] (2.1) 

where : p *= density 

pu = linear momentum in the x-diiection. 
pv = linear momentum in the y-direction. 
pe = total energy. 
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These four conservation variables can be used to calculate the following four 
primitive variables : p , u , v and p. 

where: u = velocity in the x-direction. 

V = velocity in the y-direction. 
p = pressure. 

Note that the pressure is calculated from the following formula for a perfect 

gas: 


P = (Y-l)P 


e - 


2 2 
U +V 


where : p , e , u and v are as specified in equation (2.1). 
y = ratio of specific heats for the gas/liquid 
under consideration. 


The fluid flow is governed by the balance laws of mass, momentum and 
energy conservation. These equations, expressed in integral form, arc known as the 
Euler equations for flow of a compressible inviscid fluid. This integral form can be 
written as follows : 

4-fUdn = -fQ-nds (2.2) 

an 

where : U = solution vector. 

df2, ds = Lebesgue measures of length and area, 
n = vector normal to the boundary. 
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Q = pair of flux vectors E, F. 


The form of Q is as follows : 


q = [e.f]. 


pu 

pv 

2 


pu +p 

pvu 


2 

puv 

pv +p 

(pe+p)u 

(pv+p)u 


Besides satisfying equation (2.2), U must also satisfy an initial condition 
given at t = 0 ,ie. 

U(x,0) = Uo(x) xeQ 

Now to obtain a Galerkin approximation, we must consider the differendal 
form of equation (2.2). 

^ + V. Q = 0 (2.3) 

2 9Q . 

where : V • Q = X -=r— 

dxj 

If we multiply equation (2.3) by a test function (p^ and integrate, we obtain 
the following weak boundary-value problem which is equivalent to the conservation 
laws, the initial conditions and jump conditions : 
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CHAPTERS 


FLUX^ORREOED TOANSPOirr METHO 

In Chapter 2 , a weak variational statement of the conservation laws was 
formulated, governing the flow of a fluid through a re^on Q (see equation 2.4). For 
such flows, point and line discontinuities can occur in the' primitive variables in use. 
Depending upon the scheme being used, these discontinuities can cause problems. 
In particular it is known that those schemes of order greater than one will tend to 
cause oscillations in the solution at and about such discontinuities. If these 
oscillations are large enough, the solution will eventually become unstable and 
diverge. A method able to deal with these osdillations would hence enhance the 
solution greatly. 

t 

The Hux-Corrccted Transport (FCT) Method is such a method as it employs 
the use of both a high- and low-order munerical scheme. The idea behind this is to 
use the high-order scheme in areas where the primitive variables change smoothly 
and not abruptly. The low-order scheme is then employed in those areas where the 
variables vary abruptly (such as along the line of a shock-wave). The combination of 
these two sqhemes near such discontinuities tends to provide an accurate picture 
although slight oscillations can still mar the solution. For this reason, a strong 
diffusion term is added to the low-order solution which tends to reduce the 
magnitude of these spurious solution oscillations. 
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Before a detailed description of the high- and low- order schemes is given, a 
rough outline of the FCT method is given: 

If we discretize equation (2.3) with respect to time, we obtain an equation of 
the following form (using the standard high-order method) : 

U'^^ = U" + AU*' (3.1) 

where : A s the increment in the solution vector corresponding 

to a change in time firom to tn^j. 

The FCT method computes a A U** of high enough order to capture the 
solution with few oscillations. Rewriting equation (3.1) in terms of low- and high- 
order contributions yields the following : 

_,n+l _.n . _,1 . . _,h . _, 1 . 

U =U +AU +(AU -AU) 

= + (AU** - AUS 

where : U* = low-order solution at t=tn+i. 

A U* = low-order increment in U from t=t„ to t=tn+j 


Hence the FCT method limits the second term above and sets it equal to AU* 
,,n+l ■ ,1 ^ • 

U = U + A U (3.2) 

where : A U* = ( A U*’ - A U* ) 
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Obviously A U* must be "limited" very carefully so as to avoid either 
oscillations or lack of resolution in the solution. Before describing any finer details 
of the FCT method, the high- and low- order schemes will be described. 


3.1 High-order Scheme : Tavlor-Galerkin method. 


First we discretize U with respea to time using the midpoint formula : 


,,n+l _,n _,At., 3U . 3^ 

U = U + — ) + 0(At ) 


_,n ^ 

= U + At(-^) 


(3.3) 


where : U = solution vector at t=t"+l. 

U " = solution vector at t=t " . 

( 9U / 9t solution derivative at t=t"+^/2_ 

From Chapter 2 and the governing differential equation (equaton 2.3) : 


-I' 


(3.4) 


( obtained by substituting :V-Q = V- [E ,F] 


^ 9F \ 
9x ^ 9y ' 
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If equation (3.4) is substituted into equation (3.3), the variational form of 
equation (3.3) becomes: 



Now, we know that : 


n+— X n+-i 

^,E = E = or 




’!i > 


1 

rH— ocp; 

I* 2 J 


dx 


Similarly fwF ; 


( T 


n+— *p 


a -- 2 - 

>i) 


I>f 


± dq>, 

2 T) 
dy 
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Using the above, equation (3.5) becomes : 


f J f ¥1^ ^ J 

J U dx = J U dx 


n T 


nf j 3<Pj nfy 


r 3 T 3 *^7 T 

_AtJ(^(E ^9j) + ^(F ^9j))dx 


"i 


+ AtJ(E 

a 


nf| 3q>: nf- 3q); 

+ F ^V-)dx 
dx dy 


-At](E ^n + F ^ n^ )ds 

an ^ 

where : n^^ , ny = rectangular vector components of the 

normal to the boundary (see Figure 2.1). 
an = boundary of Cl (see Figure 2.1). 

Finally, the solution vector U is discretized by : 


(3.6) 


, NNOCe , 
T,n+1 Y m-1 T 
U = L Uj q., 


NNODE 


,,n Y 

U = 2. U; <Pj 


t=l 


where : NNODE= number of nodal points. 
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9i = Galerkdn shape functions (provide 
a basis of finite dimension). 

Equation (3.6) thus becomes : 


NNODE . _ _ NNODE^ _ _ 

§ “i <Pi 15 = 2, Ui ^ <Pi 15 <h 


, , ,, 


Bx 


f !>+•*“ T T 

-AtJ(E n^tpj + F ny<pj )dx 


X X 

If we denote^ dx as Mj j , the mass matrix, equation (3.6) can be written 


as follows : 


M M U" + Q 


(3.7) 


iM-i-3q>: nf^ 3q>: 

where: Q = At|(E + F " 3 ^)d^ 


1 

-AtJ(E ^n^q)j + F ^Hyq)j )ds 
dCl 


Note that Q is expressed in terms of and F"'''!/^ which are both 

functions of .Everything else in equation (3.7) is known or can be computed 
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(c.g. ) if can be computed. To calculate , we again use equation 

( 2 . 3 ), formed in Chapter 2. The process is identical to the one for computing 
except for the time interval which spreads from t=t” to t=t’”‘*^. Hence we use : 


= U" + y ( ^ ) + 0(At^ 


u" - — ( e" + f" ) 


The following weak statement is made: 


r>f — 


j u" ' (p/ dx = (p/ dx - + Fy)(p/ dx (3.8) 


At f.^n „n, T 


U" (and hence E;j , Fy ) are known, enabling to be computed. A 

brief summary of the high-order process is as follows : 

(1) Using equation (3.8), compute knowing U", E" and F". 

(2) Substitute into equation (3.7) and compute 

(3) Repeat process to advance a further time-step. 

3.2 Low-order scheme : Lumped mass matrix. 

As mentioned earlier, the low-order scheme is used really to capture line or 
point discontinuities. As such, it should provide an oscillation-free solution to 
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prevent possible later numerical instability. The scheme uses a lumped mass-matrix 
Ml instead of the matrix used in the high-order scheme (M)- Mi is obtained by 

s ummin g the elements of each row and placing the sum in the diagonal position ,e.g. 


'5 

0 

0 

0 

O' 


'4 

1 

0 

0 

O' 

0 

4 

0 

0 

0 


1 

2 

1 

0 

0 

0 

0 

3 

0 

0 

M = 

0 

1 

1 

1 

0 

0 

0 

0 

4 

0 


0 

0 

0 

2 

2 

0 

0 

0 

0 

5 


0 

0 

0 

3 

2 


The form of the resulting equation is similar to equation (3.7) except for the 
addition of a strong diffusion term. This diffusion term, V , aids in dampening out 
any oscillations in the solution near the discontinuities (shocks) and is included as 
follows : 

jykl 

M^Uj =MiUj-hQ-»-V (3.9) 


3.3 Iterative form. 


The low- and high- order schemes have the following form : 

Q (high-order) 

MjU =M|U + Q + V (low-order) 
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Equation (3.7) can be rewritten as : 


MU +^U -MiU = MU + Q 

n+1 n4>l n 

A (M-Mj)U + MiU = MU + Q 

The iteration is as follows : 

Ml U^""^ = MU" + Q - (M - Mi)U^,5 
where: U^^ = "i ’th" iteration of 


(3.10) 
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3.4 Flux limiring procedure ; 


If we take equations (3.9) and (3.10) and subtract them, we get : 

- uJ^S = (M-l^)u"-(M-Mi)U^,S - V or 

+ (M-Mi)U" 

1x4.1 

-(M-Mi)Ujr,j - V or 


U 


n+1 

0 


U 


n+1 



which is similar in form to equation (3.2). The limited increment A U* then 
corresponds to fy.i]. 

3.5 Error estimate. 


As the "solver" based on the FCT-FEM method was to be used with 
adaptivity (see Chapter 4), an error estimator, capable of identifying abrupt changes 
in the solution vector U, was needed. The error estimator had to be accurate enough 
to control refmement/unrefinement in areas in the mesh of large/small error. Many 
references were available in the area of error estimates (see [1 1], [12]) enabling an 
appropriate error estimator to be chosen. The form of this estimator is given by the 
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normalized gradient of one of the four conservation variables. In this case density 
was chosen : 




6e = 


max 

i»U 




Ph 


where : = area of the element 

Pjj = average value of the density over element Cl^ 


3.6 Implementation of the Flux-Correc ted Transport method. 


An existing FCT-FEM "solver" using quadrilateral elements was used (see 
[13]). As this "solver" was written with quadrilateral elements in mind, it had to be 
converted in order to use triangular elements. This involved going through the code 
and changing all loops and subroutines involving elemental calculations (e.g. number 
of sides/nodes per clement changed from 4 to 3). The common blocks and 
variables/arrays were then used as the basis for the adaptive part of the code (see 
Chapter 4). This ensured a convenient way of linking the FCT-FEM "solver" and the 
adaptive algorithm described in Chapter 6. 
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CHAPTER 4 


MESH REFTNEMEm* STRATEGY 

Many reflnement strategies presently involve meshes consisting of 
quadrilateral elements. These strategies have to use constrained nodes (ie nodes 
generated by the refinement whose nodal values arc constrained by the nodal values 
of the two nodes .on the same element side) . These constraints have to be 
incorporated when computing a finite element soludon and hence result in a more 
complicated and less efficient code in general. When considering a triangular finite 
element mesh, however, one can devise refinement methods involving two 
alternative element divisions which would eliminate the need for constrained nodes. 
Such a method would thus improve the speed with which the code computed the 
solution. This is important as the majority of the running time is spent on computing 
a solution rather than on refining/unrefining the mesh after a specified time interval. 
This indicates that the more efficient and quicker adaptive codes involve triangular 
elements and unconstrained nodes. 

Taking the above into account, an adaptive triangular element strategy was 
adopted. Two different clement sub-divisions arc used in this strategy so certain 
unsatisfactory refinemcnt/unrefinemcnt arrangements can occur unless panicular 
rules arc devised and implemented. For the sake of case of reading, the type of sub- 
divisions will be described before the rules governing refinement/unrefinement. 
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4.1 Tvp>es and uses of each division type 


The two types of division, following the example described in [1], are 
known as "regular" division and "green" division. The "regular" division results in 
the generation of four identical triangular elements which are all geometrically similar 
to the original "father" element. This node is then connected to the vertex on the 
opposite side, thus dividing the element into two unsimilar elements.(See Hgure 4. 1 
for a schematic representation) 

The way in which these alternative divisions are used is as follows: 

(1) Based upon some a-posteriori error estimate, "regular" division is carried 
out throughout the mesh in such a way that there are no other elements with 
more than one constrained node per side, i.e. after "regular" refinement, the 
mesh consists of only triangular elements and degenerate quadrilateral 
elements (degenerate because of the fourth node introduced by refinement- 
nodes 1-3 are such nodes for the shaded elements in Figure 4.2 a ). 

(2) "Green" refinement is then carried out throughout the mesh on all of the 
degenerate quadrilateral elements. This effectively "cleans up" the mesh and 
yields no elements with constrained nodes (refer to Figure 4.2 b ). 
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NODES 



a) "Regular” division 


NODES NODES NODES 





b) "Green" division 


Figure 4.1 : Types of division 
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4,2 General Rules governing Refinement/ Unrefineinent 


As the refinement/unrefinement must produce reasonably shaped elements 
without constrained nodes, two general rules can be immediately defined: 

(1) The interior angle of each element must be bounded from zero to ensure a 
reasonable soludcm 4>e. long, slender elements are not desirable. 

(2) The size difference between neighboring elements must not be such that 
constrained nodes are produced. 

4,3 Rules governing successive refinements 

Rule number (1) above prohibited long, slender elements. The combination 
of bodi "green" and "regular" divisions, however, can contradict this rule as can be 
seen in Figure 4.3. These figures show situations which could contradict rule (1) 
and the solutions to these situations. 

Rule (2) dissallowed constrained nodes. Figure 4.4 shows situations which 
could lead to at least two constrained nodes if care is not taken in refinement. 
Depending on the element numbering, the "cleaning up" of the mesh (see Figure 
4.2) with "green" divisions would reduce these two nodes to one constrained node 
which is still unsatisfactory .The following logic is used to solve these problems: 
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o : constrained nodes resulting from 
incorrect refinement 



(b) 




Figure 4. 4 Funher examples of correct and incorrect successive refinements 
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NOTE: In the following, NELl=the element that is to be refined due to 

a large a-priori error estimate. 

NEIG= any neighbouring elements of NELl 
NNEI= any neighbouring elements of NEIG 
LEVEL(NEL)= level of refinement of NEL 

( =0 for a non-refined element) 

(1) If NELl is a "green" triangle , NEL must be "un-greened" and 
then"regularly" divided before NEL can be divided( Figure 4.3 a ) 

(2) If any NEIG is a "green" Triangle, NEIG must be " un-greened" and 
"regularly" divided bfore NELl can be divided( Figure 4.3 b ) 

(3) If any NEIG has one or more "regularly" divided neighboursQ^NEIG), 
then NEIG must be "regularly" refined before NEL can be dealt with (see 
Figure 4.4 a ). 

(4) If LEVEL(NEL) is greater than LEVEL(NEIG), NEIG must be refined 
"regularly". (Figure 4.4 b ) 
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4.4 Rules eovcming successive unrefinements. 

In the case of unrefincment, the two general rules (see section 4.2) must 
again be obeyed. If refinement has been carried out in the correct manner, 
unrefinement cannot cause long, slender elements (rule (1)). Hence rule (1) will 
automatically be obeyed in unrefincment. Constrained nodes (rule (2)), however, 
can occur if certain rules arc not obeyed. One clear case can be seen in Figure 4.5 
where removal of the shaded group would cause two constrained nodes to occur, 
destroying the continuity of the solution. The following logic is hence used: 

NOTE: In the following: NG1= group number requiring unrefincment. 

NEIGRP=any neighbouring group of NGl. 

NNEIGRP=any neighbouring group of NEIGRP. | 
"REGULAR" group : group of "regular" elements. 
"GREEN" group : group of "green" elements. 

t 

(1) If NGl is a "green" group, unrefinement is not carried out. 

(2) If NGl contains a further group (NG2), group NG2 must be unrefined 

before NGl is considered again (see Figure 4.6) . 

(3) If two or more of the three NEIGRP's are "green" groups, unrefinement 

is immediately possible (see Figure 4.7). 
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(4) If two or more of the three NEIGRP's are "regular" groups, then the 
three -NNEIGRP's of these NEIGRP's must be examined: 

a) If two or more NEIGRFs have one or more "green" NNEIGRP's 
uiuefinement is carried out (see Figure 4.8). 

b) If one or more NEIGRP's have one or more "green" groups as 
neighbours and one NEIGRP is a "green" group, unrefinement is 
possible (see Figure 4.9). 

c) If neither situation a) nor b) exists, then unrefinement is deemed 
impossible. 

4,5 Pnctical implementation of Refinement/Unrefinement 

Before describing the subroutines which control refinement and 
unrefincment, a brief outline of the data structures used is necessary. The data arrays 
necessary for refinement/unrefinement will be explained in detail as well as the 
original data structures contained within the finite element "solver". 

4.5. 1. Data structures used: 

The original code (ie the "solver") needs the following data structures to 
provide a solution: 

NELEM= number of elements. 

NNODE= number of nodes. 

NODES(I,NEL)= Global node number of the I’th node of element "NEL". 
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urine "regular” groups C hav 







X(IJNOD)= I'th coordinate of global node number "INOD". 

QTN(I4NOD)= rth component of the solution for global node "INOD". 

The following arrays were introduced for refinement 

NELCONflJ^EL)® Fth connection/ncighbor of element "NEL". 
LEVEL(NEL)= level of refinement of element "NEL". 

ERR(NEL)= an a-postcriori estimate of the local error of the solution in 
element "NEL". 

A brief explanation of the NELCON and LEVEL arrays follows: 

(1) NELCON array : this array stores the neighbours of a particular element 
The connectivity numbering of an arbitrary element is shown in Rgure 4. 10. 
The dimension of NELCON is (6.MAXEL) as each side of the element can 
be connected to two different elements( MAXEL = maximum number of 
elements ). 

(2) LEVEL array : every time an element is refined, the level of the resulting 
two or four elements ("green" or "regular" division) is increased by one. 

For unrefinement to occurr, we need three additional data arrays . One of 
these is used to store some history of the refinement process and the other as an 
error indicator for each group of elements generated by refinement (unlike 
ERR(NEL) which is specified per element). 

NELGRP(I,NG)= I'th member of group number NG. 

If > 0 : Refers to an element 

If < 0 : Refers to a group of elements. 
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GRERR(NG)= an error estimate for groups of elements generated by 

refinement (GRERR= the sum of ERR's of each element in 
the group). 

IELGRP(NEL)= the number of the group containing element NEL. 

For a schematic representation of the NELGRP array, see Figure 4.11, 
noting the differences between "green" and "regular" groups. Lastly, in order to 
control the generation of new nodes, elements and groups, an additional structure 
obtained from [2] was implemented. 

INODFR(I)= an array containing a list of "free" nodes, ie nodes not 
currently being used, "freed" by uruefinement . 

INELFR(I)= an array containing a list of "free" elements caused by 

unrefinement (as opposed to INODFR(I) which contains "free" 
nodes) . 

IGRFR(I)= an array containing a list of "free" groups caused by 
unrefinement 

IELCH(I)= an array containing a list of elements which have either been 

refined or unrefined and are also still in use (as opposed to "free" 
elements). 
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ELI 


GROUP NGl 



NELGRP(1,NG1) = ELI 
NELGRP(2,NG1) = EL2 
NELGRP(3,NG1) = EL3 
NELGRP(4,NG1) = EL4 


a) "Regular" group 



NELGRP(2,NG1) =.EL1 
NELGRP(3,NG1) = EL2 
NELGRP(4,NG1) = 0 
b) "Green" group 

Fieure 4. 1 1 Schematic representation of NELGRP array. 
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4.5.2 Subrourines controlling Refinement/Unrefinement. 


4.5.2. 1 Refinement. 

Refinement is obtained mainly through the use of subroutines REFINE and 
DIVIDE(NEL1,NEL2). REFINE decides which elements need to be refined based 
on our error array ERR(NEL) and provides DIVIDE with these element numbers 
(see Figure 4.12 for a flowchart). Subroutine DIVIDE then operates as follows: 

SUBROUTINE DIVIDE(NEL1.NEL2) 

NELl : (input) 

NEL2 : (output) =NEL1 if NELlhas been been refined. 

= NELO if NELO has to be refined first to enable 

NELl to be refmed.All the rules outlined 
in section 4.3 are implemented in DIVIDE. 

Once all the elements have been considered for "regular" refinement, a final 
loop over the elements (in REFINE) identifies any degenerate quadrilateral elements 
and "green's" these elements to ensure no constrained nodes. Note that once an 
element is to be divided, the following steps are necessary: 

(1) Generate new element numbers. 

(2) Change NELCON array of neighbours. 

(3) Generate NELCON array for new elements. 
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FOR EL=1,NELEM : 
IS ERR(IEL)>EMAX? 


NEL1=IEL 



o^p : NEL2 


SUBROUTINE 
DIVIDE (NEL1.NEL2) 


IS NEL2=NEL1 ? 



IS 1=0? 


1 = 1+1 

NEL1=NEL2 

NLa)=NELl 


NELl=NLa) 

1 = 1-1 


Figure 4.12 Flowchart of how subroutine REFINE calls subroutine DIVIDE 
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(4) Generate new nodal numbers. 

(5) RU in the NODES array for the new elements. 

(6) Interpolate the soludon over the new nodes. 

4.S.2.2 Unrefinement. 

UnreHnement is obtained through the use of subroutines REFINE, 
UNREFINE(NG1,NG2) and (sec Figure 4.13). REFINE again 

decides which groups require unrefinement (based this dme on the group error array 
GRERR(NG)) and provides UNREFINE with these group numbers as follows. 
Subroutine UNREFINE then operates as follows: 

SUBROUTINE UNREFINE(NG1,NG2) 

NGl : (input) : group to be unrefined. 

NG2 : (output) = NGl if NGl has been successfully unrefined. 

= NGO if Group NGO has to be unrefined first in 
order to unrefine NGl (see Section 4.4 for 
unrefinement rules). 

= 0 if unrefincment of group NG 1 is impossible. 

NOTE : that if NG2 has a group error which does not demand unrefinement, then 
NGl will not be unrefined although its group error does require unrefinement. This 
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is unlike Subroutine DIVIDE(NEL1,NEL2) which will refine any elements just to 
divide NELL 

Again once a group is to be unrefined, the following steps are required: 

(1) Select one element number from the group (of two to four elements). 

(2) Adapt the NODES array of the neighbouring elements. 

(3) Adapt the NODES and NELCON arrays for the selected element. 

(4) Add one or three elements to the list of free elements depending upon 
whether the group was "green" or "regular". 

(5) Add the group number to the list of free nodes. 
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CHAPTERS 


NUMERICAL EXAMPLES 

In this section, examples of both steady-state and transient problems -in 
supersonic compressible flow are presented. The Flux-Controlled Transport (FCT) 
algorithm together with the adaptive refinementAinrcfincmcnt algorithm was seen to 
produce very good results especially in those problems involving shock waves (i.e. 
line discontinuities) and Pfandtl-Meyer expansion waves. In almost all of the 
examples exact theoretical solutions were available to compare with the computed 
solutions, allowing direct comparisons. The FCT algorithm was employed in all the 
problems beside problem 5.2 with a value of c = 0. 125. 

5.1 Supersonic flow over a 20* concave comer (’ramp). 

In this case, supersonic fluid flow over a 20’ ramp was observed. When 
supersonic flow is deflected_upwards through an angle 0 , the flow streamlines have 
to change direction very abruptly. This takes place across the shock wave which is 
oblique to the initial flow direction and stems from the point at which the flow is 
deflected. All the deflections are alike meaning that the flow remains parallel after the 
shock. Across the shock, the fluid velocity decreases and the density, pressure and 
temperature all increase. Refer to Figure 5.1 for a sketch together with the initial 
coarse mesh. 
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Two levels of refinement were allowed and the constants controlling 
refinement/unrefinement were p = 0.20 and a = 0.01, i.e. the error estimate for each 
element ideally lies between a*EMAX and P*EMAX (EMAX is the maximum 
element error estimate, selected by looping over all the elements) . Uniform inflow 
conditions of Mach 3 and 7=1.40 ( Y=Cp/Cv) were specified. After less than 500 

time steps, the shock profile had been accurately captured except for some small 
disturbances which ^sapeared ^ter further time-stepping. The solution at this time 
is shown in Figure 5.2 together with the final refined mesh. The concentrated area 
of level 2 elements (level 0 = unrefined elements) correlates excellently with the line 
of the shock. The slight oscillatory nature of the solution is characteristic of the FCT 
method. 

5.2 Supersonic flow over a convex 10* comer. 

In contrast to flow over a concave comer (section 5.1), flow over a convex 
comer results in the fluid being deflected away from itself to remain parallel with the 
surface. This change in direction is accomplished through an expansion wave which 
is centered at the point at which the comer begins. The flow streamlines are all 
uniformly curved by the expansion fan until they are again parallel to the surface. 
Unlike the discontinuities across a shock wave, the flow properties change smoothly 
and continuously over an expansion wave. In addition the flow velocity increases 
whereas density, pressure and temperature all decrease (see Figure 5.3 for a sketch) 

Two levels of refinement were allowed together with the following 
refinement/unrefinement constants 7=1.38 , P=0.20 and a=0.01. The inflow 
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(a) Rnal computed solution 



(b) Final refined mesh (maximum level of refinement =2) 
Figure 5.2 Final soltirion and the corresponding refined me<;h 
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conditions were specified as a uniform Mach 6 flow. Various values of c were 
allowed ranging from c= 0.0 to 1.0 . The difference in the amount of diffusion 
correlates to the value of c as can be observed from Figure 5.5. The final solution 
(c=0.125) and the corresponding refined mesh are shown in Figure 5.4. 

Again the results compared very well with the exact solution. The area of the 
expansion fan is made up from level 2 elements as expected. Because of the high- 
speed inflow conditions, a shock wave is generated at the second (concave) comer 
after the expansion wave. This is as expected remembering that the fluid velocity 
increases through an expansion wave, ensuring continued supersonic (if not 
hypersonic) flow. 

5.3 Intersection of two shock waves of the same family. 

A sketch of the double-ramp problem (appropriate to this class of problems) 
is given in Figure 5.6 as is the initial coarse mesh. The supersonic inflow conditions 
are specified so as to generate two different oblique shock waves (similar to problem 
5.1) , one at A and one at B. Shock wave BC, because of the increased ramp angle 
at point B, will be inclined at a steeper angle than shock wave AC. Hence the two 
shock waves intersect at point C resulting in the propagation of a single shock CD. 

Now the flow direction and pressure in region 3, P3 and 63, result from the 
upstream conditions in area 2.Likewise p2 and 02 are a result of upstream conditions 
in area 1, so P3 and 63 are affected by both shocks AC and BC. Propenies in area 5 
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Figure 5.4 Final solution and the cnrre.sponding mrch 







•weak reflected 


(a) Intersection of two shocks generated by a double ramp. 


mjd 


(b) IniuaJ coarse mesh 

of two shocks intersecting and the initial mesh 
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(P5 and O5), however, are processed by a single shock CD. As the entropy changes 
across the single shock (CD) and the two shocks (AC and BC) will be different, a 
slip line must originate from point C (ie a line across which the pressures and flow 
directions are equal). To get P5=P3 and 05=63 simultaneously is virtuaUy impossible 

and hence a weak reflected wave is generated from point C (either a weak shock or 
an expansion wave). All this wave really does is to ensure that P4=P5 and 64=65, 

satisfying all shock relations. 

Three levels of refinement were specified with y=1.38, P= 0.20 and a=0.01 . 
The inflow conditions were uniform Mach 5 flow with c=0.125. 

The results as can be seen in Figure 5.7 correspond excellently to the 
expected appearance (that contained in Figure 5.6). Two shock waves are generated 
and intersect, combining to form a third wave. A weak shock wave can be seen to 
originate from point C to satisfy all physical relations between area 1 and area 5, and 
between areas 1, 2, 3 and 4. The area of most refinement is concentrated about all 
four shock waves (three strong waves and one weak wave) as was expected (see 
Figure 5.6). 

5.4 Shock Reflection problem. 

This problem is generated by specifying initial conditions corresponding to 
four different triangular areas. It is imponant that these initial values of density, 
momentum, and energy satisfy the Rankine-Hugoniot jump conditions over a shock 
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wave. Although this problem might seem a little contrived, the solution is identical 
to, say, that of a conventional supersonic diffuser. In this case, the domain is a lot 
less complicated (see Figure 5.8). 

The inflow conditions correspond to a Mach number of 3.5 and y= 1.4. The 
refmement/unrefinement parameters are again (3=0.20 and a=0.01. 

Figure 5.9 is a plot of the interpolated initial conditions and the mesh that 
results from initial refinement (before time-stepping has begun). The computed 
solution (after 500 time-steps) is then shown with a further refined mesh in Figure 
5.10, One can observe that more elements have resulted from additional refinement. 
The mesh compares favourably with those obtained in [12] as does the converged 
solution. The FCT algorithm again gives the solution an oscillatory appearance 
which does not detract from the overall appearance of the solution. 

5.5 Supersonic flow over a step. 

This problem is an^agous to flow in a wind tunnel or a long tube over a 
small step (in this case the step height was a quarter of the tube diameter). The mesh 
and solution plots have all been intentionally scaled incorrectly for reasons of detail. 
The initial mesh in both scales is given in Figure 5.11 to illustrate this. 

The inflow conditions were specified as Mach 3 and y= 1.4 on the left venical 
mesh edge (the right vertical edge was specified as outflow). All other edges 
correspond to no-flow conditions. 





(b) Reflected waves inside a diffuser 

Figure 5.8 Initial conditions for reflected wave problem and a physical example. 







(a) Initial interpolated conditions 



; ^ (b) Initial refined mesh (time = 0 sec. ) 

Hgure 5.9 Initial conditions and mesh for reflected wave problem 

J 


1 


J 


D-53 













HKBbbbbbSi^— - 




Mi^SJJSSSSSm’MSSMiJi! 


HHaBBggSSiiiigggMBHBI 


BBBBBBS 






SBBBBBBB 


HBHBBHBiiiiiiSBHBBHHB 



WlMrArA^^^ 




mw&ws 


MBIHBIHHBIBBBBBF^K8gs»«i 





(b) Final refined mesh (time = tfinal) 

Figure 5,10 Final solution and mesh for reflected w.ave problem, 
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The solution to this problem involves several shocks/ shock interactions and 
took 1000 time-steps to converge. The problem was run for two cases of different 
MAXLEV (= maximum level of refinement for an element). The final converged 
solutions and the corresponding refined meshes are included in Figures 5.12 and 
5.13. 


Figure 5.12 (MAXLEV=2) shows a similar shape solution to that contained 
in Figure 5.13 (MAXLEV=3) although the shock’s resolution is better defined and 
clearer in Figure 5.13. The basic form of the solution involves an upstream bow 
wave caused by the step. This wave then reflects off first the upper surface then the 
lower before exitting. A Prandtl-Meyer expansion wave is generated at the comer of 
the step. The areas of large error (corresponding to the areas contining the shocks) 
arc again well "covered" by the smaller (more refined) elements. Because of the 
complexity of the solution, the oscillatory nature of the FCT method is clearer 
especially in the second and third (reflected) waves. 
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(a) Final computed solution 



(b) Final refined mesh 
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(b) Final refined mesh 

Bgurc 5.13 Final solurion and correspondine refined mesh rMAXLEV=3) 
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CHAPTER 6 


RECOMMENDATIONS AND CONCLUSIONS 

The main aim of this report was to provide a code capable of providing high- 
resolution solutions to a variety of problems involving strongly unsteady 
compressible flow. This was to be done with the minimum number of elements 
necessary. 

As far as these two objectives are concerned the code performed very well, 
providing both accurate pictures of shocks and shock interaction as well as optimal 
meshes (due to the refinement/unreflnement capabilities of the code). The 
combination of the FEM-FCT numerical scheme and adaptivity was seen to be well 
suited to the class of problems of analysis. All the solutions were true to the 
calculated or experimentally-obtained exact solutions. 

The success of this scheme would seem to prompt the development of a 
three-dimensional code. There are several references to the multi-dimensionality of 
the Flux-Controlled Transport method (see References [14], [10]). The types of 
refinement/unrefinement could also be generalized to three-dimensions. In the case 
of "green" divisions (see Chapter 4), the three-dimensional equivalent would 
generate four tetrahedral elements. Similarly, "regular" division would cause the 
generation of eight geometrically similar (to the parent element) tetrahedra. The 
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combination of these two divisions would result in three-dimensional meshes with 
no constrained nodes, again indicating possible time savings in the solution 
computation. Three-dimensional solutions to the integral form of the Euler equations 
would have great use and potential, providing information about the effects of 
surface shocks in three dimensions. 
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Appendix E 


This appendix documents the implementation of an adaptive strategy using 

Lockheed-Huntsville GIM/PAGE code as the flow solver. A concise and complete 

discussion of GIM/PAGE methodology and its application to practical flow 

* 

problems can be found in the work of Spradley et al. 

The implementation of the adaptive strategy using the GIM/PAGE code is 
possible because the adaptive strategy does not have to be used with any 
particular solution algorithm. For this reason, many existing CFD codes can 
be made adaptive without a major effort. The CFD code must be able to properly 
treat •'constrained” nodes which will exist along the interface between the 
different levels of refinement. If this capability is not present it must be 
added. Routines must be added to transfer information between the CFD code 
and the adaptive strategy. 

A general adaptive strategy for the computation of steady-state solutions 
of the equations of compressible gas dynamics involves the following steps: 

1. For a given mesh, compute the steady-state solution. 

2. Compute the local error for the mesh. 

3. Survey the error field and determine where mesh restructuring 
is needed. 

4. Refine or coarsen the mesh as needed. 

5. Go to step 1. 

The general tendency is to combine the adaptive strategy and flow solver 
into one large, complex computer program. However, a closer examination of 
the general adaptive procedure suggest that the development of another large 


^Spradley, L.W., Stalnaker, J.F., Robinson, M.A., and Xiques, K.E., ’•Finite 
Element Algorithms for Compressible Flow Computation on a Supercomputer," 
Finite Elements in Fluids (eds R.H. Gallagher, G. Carey, J.T. Oden and O.G. 
Zienkiewicz) , Vol. 6, 1985. 
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computer program can be avoided. The Lockheed-Huntsville PAGE code can already 
perform step 1 of the adaptive strategy. Furthermore, steps 2 through 4 do not 
require the resources of a supercomputer. They could be executed using the 
less expensive resources of a front-end machine or workstation. 

This adaptive strategy using the GIM/PAGE code has been implemented 
through a Solution-Adaptive Analysis System (SAAS) . This system operates in a 
computational environment consisting of a supercomputer and a front-end 
machine. The SAAS concept is shown in Fig. E-1. The adaptive processor and 
the adaptive database together with pre- and post-processors used to transfer 
data to and from the PAGE code reside on the front-end side of the SAAS 
environment. The PAGE code is used on the supercomputer side of SAAS. This 
segregation of the mesh restructuring from the flowfield computation will 
allow SAAS to be used with other existing CFD codes. Only custom pre* and 
post processors need to be developed to facilitate data transfer between the 
adaptive processor and the CFD codes. Currently, only pre- and post- 
processors which interface with the PAGE code are available. 

The SAAS performs solution-adaptation through its adaptive processor. 
Utilizing the current geometry and solution, the adaptive processor first 
assesses the solution quality over the entire domain. It then refines or 
coarsens the mesh as required based on user-supplied refinement and coarsening 
tolerances . When refinement (coarsening) is indicated, the physical grid 
spacing in each of the computational coordinate directions is halved (doubled), 
the local nodal connectivity is reconfigured, and all geometric and flowfield 
variables are updated. Flow field quantities at grid points which fall on the 
boundaries between refined and non-refined elements are constrained to values 
determined by linear interpolation between the associated connected nodes. 

The pre- and post-processors convert the geometry and solution data to and 
from PAGE code format. The PAGE code is then used to compute a new solution 
using the adapted geometry. The SAAS work environment creates no computer 
resource penalty during the actual flowfield integration. This is because the 
grid refinement is removed from the integration process and is performed on a 
less expensive workstation or front-end computer. 
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The following test cases were run to determine the utility of the SAAS 
adaptive work environment. The refinement threshold was held at 0.65. First, 
the flow of an ideal Mach 2.4 freestream up a 14.04 deg ramp was computed on a 
relatively coarse grid using the PAGE code. The grid and resulting solution 
are shown in Fig. E-2. SAAS was used to adapt the grid and compute an en- 
hanced solution using the PAGE code. The adapted grid and the subsequent PAGE 
code solution are shown in Fig. E-3. 

Second, a viscous shock-expansion test case was run on the same config- 
uration and with the same inflow. A Reynolds number of 1000 was used. SAAS 
was used to adapt the mesh four times. The PAGE code was used to compute an 
intermediate solution between each adaptation. Figure E-4 shows the final 
grid and subsequent PAGE code solution. 

Figures 5 and 6 show an inviscid shock reflection calculation. The con- 
ditions were the same as the previous cases with a solid upper boundary. This 
case was run to test the ability of SAAS to capture and refine multiple inter- 
actions of shocks and expansions. The mesh shown in Fig. E-5 was refined 
twice yielding the pressure results shown in Fig. E-6 . 
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Appendix F 


The three-dimensional adaptive scheme is divided into subprocesses each 
of which performs a specific task or function. The subprocesses are: 

1. Feature detection 

2. Mesh adaption 

3 . Boundary condition application 

4. Data translation 

These subprocesses surround a central data pool as shown in Fig. F-1. This 
arrangement is essentially the same organization used by Dannenhoffer and 
Baron (Ref. 1) in their hybrid expert system. 

Two-way data transfer is allowed between each subprocess and the central 
data pool. There is no communications between any two subprocesses. As 
Dannenhoffer and others have previously stated, this type of communication 
structure benefits the overall system design and development. Each subprocess 
may be developed and tested separately. Each subprocess is independent of the 
other subprocess . Any algorithm change in a subprocess will not impact any 
other subprocess. All requests for data by any subprocess is handled by a 
collection of procedures which read and write data to the data pool. The data 
pool consists of a collection of files which can reside on more than one 
directory or machine. 

The data manager processes information transfer requests between any sub- 
process and the data pool. The data manager keeps track of available elements 
and groups, retrieves information about a particular element or group and 
checks for non-existing data request. 

The feature detection subprocess determines which regions of the mesh 
should be refined or unrefined. The mesh adaptation subprocess does the 
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mesh refinement and unrefinement. The boundary condition application sub- 
process applies boundary conditions on the adapted mesh. The data translation 
subprocess converts the data to a form which can be used by a CFD code. 

The rest of this appendix describes the mesh adaptation subprocess. The 
other subprocesses have not been completed except for the data manager. 

Mesh Adaptation 

The mesh adaptation subprocess refines and unrefines the mesh using a list 
of instructions generated by the feature detection subprocess. This list 
resides in the data pool. 

Two pieces of information are necessary to adapt the mesh. These are the 
relationship of an element with its neighbors (connectivity array) and the 
lineage of the element (group array). Using the connectivity and group arrays, 
the mesh can be adapted. 

Connectivity Array 

A connectivity array is employed within the integer space to facilitate 
mesh refinement and unrefinement. This array keeps track of the dynamic 
topological relationships between the elements of the mesh. Mesh adaptation 
is performed by operating on this connectivity array. 

The three-dimensional connectivity array consists of a 4 x 4 x 4 array. 
Each element of this array is referred to as an atom. The array consists of a 
group of core atoms surrounded by neighboring atoms and neutral atoms. A con- 
nectivity array is illustrated in Fig. F~2. There is one connectivity array 
for each nonboundary element in the mesh. This array contains the topological 
data which relates one element with its immediate neighbors. Using this array 
the relative level of an element with its neighbor can be easily determined. 
Constrained nodes can also be identified. 
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Group Array 


The group array is used to facilitate the unrefinement process. This 
array contains the elements which are created when an existing element is 
refined. It keeps track of element lineage. 

Refinement Algorithm 


The following algorithm outlines the steps required to perform the 
refinement operation for a single element. 


Step 1. Copy the connectivity array of the element to work the array. 

Step 2. Get the next available group and initialize the new group array. 

Step 3. For each core location of the work array: 

3a. Get the next available element and place it in the core 
location. 

3b. Initialize next available element’s connectivity array. 

Step 4. Copy core locations of the work array to the new group array. 

Step 5. For each group location of the group array: 

5a. Get the element stored in that group location. 

5b. Transfer connectivity data from the work array to each 
element’s connectivity array. 

Step 6. For each neighbor location of the work array: 

6a. Get the element stored in the neighbor location. 

6b. Update the neighbor element’s connectivity array. 

Step 7. Update the old group array which the element belonged to. 

Step 8. Delete the element connectivity array and add element to the 
available element heap. 


This refinement algorithm is illustrated in Fig. F-3. A single group of eight 
elements is successively refined until it consists of 13 groups totaling 98 
elements. One can see from these figures that there is a maximum of one level 
of refinement difference between any two adjacent elements. An exploded view 
of the final grid is shown in Fig. F-4. 
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Fig. F-3a Original Mesh 
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Fig. F-3d Third Refinement Pass 




Unref inement AlRorithm 


The following algorithm outlines the steps required to perform the 
unrefinement operation for a single group. 


Step 1: Copy the group array to the core of the work array. 

Step 2: For each core location of the work array: (a) Get the element 

stored in the core location and (b) transfer connectivity data 
to the work array. 

Step 3: Get next available element. 

Step 4: Store this element in each core location of the work array. 

Step 5: For the group array; (a) Get an element for the group array, 
(b) update the connectivity array of each neighbor of this 
element, and (c) add the element to the available element heap. 

Step 6: If this is an embedded group then replace any reference to the 
group by the next available element. 

Step 7 : Add each member of the group array to the next available 
element heap. 

Step 8: Add group to the next available group heap. 


Figure F-5 is a refinement/unrefinement cycle using the above algorithms 
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